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The  theory  is  illustrated  by  two  examples.  A  simply-supported 
viscoelastic  beam  with  controllers  illustrates  the  solution  process 
in  Example  1.  The  beam  example  incorporates  the  fractional 
calculus  viscoelastic  behavior  in  a  structure  element,  while  an 
axially  deforming  rod  in  Example  2  incorporates  the  behavior  through 
a  viscoelastic  shear  force  applied  at  a  node  by  a  damping  pad.  The 
example  equations  of  motion  are  numerically  solved  using  the 
commercially-available  control  analysis  software  package  called 
MATRIX/. 
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Abstract 

The  objective  of  this  thesis  is  to  develop  a  control  law  for 
structures  incorporating  both  passive  damping  via  viscoelastic 
materials  modelled  by  a  fractional  calculus  stress-strain  law  and 
active  damping  by  applied  forces  and  torques.  To  achieve  this, 
quadratic  optimal  control  theory  is  modified  to  accomudate  systems 
with  fractional  derivatives  in  the  state  vector.  Specifically,  linear 
regulator  theory  is  modified. 

The  approach  requires  expanding  the  structure's  equation  of 
motion  into  a  fractional  state-space  system  of  order  1/n,  where  n  is 
an  integer  based  on  the  viscoelastic  damping  material  constitutive 
law.  This  approach  restricts  the  theory  to  materials  which  have 
rational,  fractional  derivatives  in  their  constitutive  laws.  An 
equivalent  first-order  system  is  then  formed  and  used  to  derive  the 
optimal  control  theory  in  linear  regulator  problems.  The  quadratic 
performance  index  used  in  linear  regulator  theory  is  used,  and 
equations  similar  to  those  derived  for  linear  regulator  problems 
using  first-order  systems  are  developed.  The  optimal  control  law  is 
asymptotically  linear  feedback  of  the  state  vector  for  time  large. 

An  equation  which  defines  the  optimal  gain  matrix  for  the  feedback 
is  derived  and  is  asymptotically  an  algebraic  Riccati  equation  for 
long  time  and  for  gain  matrices  which  are  asymptotically  constant 
for  large  time.  Since  an  algebraic  Riccati  equation  can  be  defined, 
current  solving  routines  can  determine  the  asymptotically  constant 


optimal  gain  matrix  for  large  time  for  the  fractional  state-space 
system. 

in  the  general  time  case,  no  Riccati  equation  can  be  derived 
because  of  coupling  between  the  optimal  gain  matrix  and  the  state 
vector  due  to  fractional  derivatives  of  the  control  vector.  Only  when 
time  is  large  and  gain  matrices  are  asymptotically  constant,  do  they 
uncouple  . 

The  theory  is  illustrated  by  two  examples.  A  simply-supported 
viscoelastic  beam  with  controllers  illustrates  the  solution  process 
in  Example  1.  The  beam  example  incorporates  the  fractional 
calculus  viscoelastic  behavior  in  a  structure  element,  while  an 
axially  deforming  rod  in  Example  2  incorporates  the  behavior  through 
a  viscoelastic  shear  force  applied  at  a  node  by  a  damping  pad.  The 
example  equations  of  motion  are  numerically  solved  using  the 
commercially-available  control  analysis  software  package  called 
MATRIXx. 


QUADRATIC  OPTIMAL  CONTROL  THEORY  FOR  VISCOELASTICALLY 
DAMPED  STRUCTURES  USING  A  FRACTIONAL  DERIVATIVE 
VISCOELASTICITY  MODEL 


1-  Introduction 

Control  of  large  flexible  structures  is  currently  of  interest 
because  of,  in  part,  the  possible  use  of  large  space  structures  in  the 
Strategic  Defense  Initiative.  Fine  pointing  and  position 
requirements  drive  the  need  for  fast,  accurate  control  of  a 
structure's  state.  Primary  emphasis  in  the  past  has  been  on  purely 
elastic  structures  controlled  by  various  applied  forces  and  torques 
in  feedback.  Unfortunately,  large  flexible  structures  have  large 
numbers  of  vibrational  modes  at  low  frequencies  which  are  easy  to 
excite.  The  modes  are  also  very  closely  spaced  in  frequency,  so  that 
control  of  one  often  excites  one  or  more  other  close  modes.  Thus, 
control  accuracy  is  degraded.  The  incorporation  of  viscoelastic 
materials  into  structures  holds  promise  to  improve  active  control 
performance  considerably. 

One  primary  reason  that  viscoelastic  materials  have  not  been 
readily  included  is  that  classical  models  of  viscoelastic  behavior 
are  complex.  The  classical  stress-strain  law  uses  a  sum  of  integer 
order,  ordinary  derivatives  to  describe  the  weak  frequency 
dependence  of  the  viscoelastic  behavior: 


m  -u  P  -\J 

S,  oO  _  v-1  d  £ 

bj  — =  Ee  +  2,  a. — 
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where  a  is  stress,  e  is  strain,  aj  and  b:  are  constants,  and  m  and  p  are 
integers.  More  or  less  derivatives  are  employed  to  fine  tune  the 
model  to  the  precision  required.  High  precision  requires  a  large 
number  of  derivatives,  implying  also  a  large  number  of  parameters 
(aj  and  bj).  The  large  number  of  derivatives  also  expands  the  state 
vector  which  must  be  observed  and  fed  back  for  control.  This 
resulting  calculations  with  large  matrices  are  not  condusive  to  fast 
response  systems. 

An  alternate  viscoelastic  model  uses  fractional  derivatives  of 
time  describe  the  material  response.  The  fractional  derivative  of 
g(t)  of  order  a,  where  a  is  taken  to  be  a  real  fraction,  is  defined  as 


dag(t)  =  i  d  g(t-x) 
dta  r(1-a)dtj0  x01 


,  0  <  a  <  1 


Its  Laplace  transform  is  [sG(s)-g(t=0)]/s1‘a,  where  s  is  the  Laplace 
variable  and  the  capital  letter  indicates  the  Laplace  transform  of 
g(t).  The  fractional  calculus  model  of  the  constitutive  law  is 
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For  this  thesis  b  is  chosen  to  be  equal  to  zero;  and  a  =  m/n,  where  m 
and  n  are  integers.  Such  a  model  is  commonly  called  a  "power-law" 
constitutive  model  (5).  Bagley  and  Torvik  (3:918)  have  shown  that 
this  model  is  accurate  over  four  or  more  decades  in  frequency  for 
over  130  viscoelastic  materials.  The  virtue  of  this  model  is  that  it 
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requires  only  four  parameters.  However,  it  too  requires  expanding 
the  state  vector  to  include  fractional  derivatives,  but  the  fractional 
derivatives  are  lower  order  than  those  used  in  the  classical  model. 
The  effects  of  noise  and  uncertainties  are  not  as  dramatically 
amplified  through  the  lower  order  diferentiation.  Consequently,  the 
fractional  calculus  model  appears  to  be  a  better  choice  to 
incorporate  into  a  control  law. 

The  objective  of  this  thesis  is  to  attempt  to  incorporate  the 
fractional  calculus  viscoelastic  mode!  into  optimal  control  theory, 
such  that  an  optimal  control  law  for  an  active  controller,  which  is 
optimal  in  a  sense  to  be  described,  can  account  for  and  complement 
passive  damping.  Current  state-space  optimal  control  theory  uses  a 
performance  index  to  optimize.  Indices  with  quadratic  terms  are 
commonly  used  because  energy  is  described  by  the  square  of  the 
velocity  and  displacement.  This  is  quadratic  control  theory.  A 
special  case,  linear  regulator  theory,  assumes  a  quadratic  form  for 
the  optimal  value  of  the  performance  index,  J,  i.e.  J  =  1/2  xTKx  , 

where  x  is  the  state  vector,  K  is  a  matrix  of  gains,  and  the  asterisk 
indicates  an  optimal  value.  In  this  special  case  quadratic  control 
theory  produces  optimal  control  as  linear  feedback  of  the  state 
vector.  Because  the  feedback  is  linear,  the  analysis  and  its 
implementation  is  relatively  straightforward.  Therefore,  this  work 
focuses  on  deriving  the  equations  in  the  linear  regulator  theory  for  a 
system  using  fractional  state-space. 

The  current  linear  regulator  equations  are  derived  by  the 
application  of  a  first-order  state-space  system  of  the  form: 
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(3) 


x(t)  =  A(t)  x (t)  +  B(t)  u(t) 
where 

x(t)  ■  state  vector 

x(t)  =  first  derivative  of  the  state  vector  with  respect  to  time 

u(  t)  =  vector  of  control  forces  and  torques 

A(t),  B(t)  =  matrices  of  functions 
The  approach  taken  in  this  thesis  is  to  formulate  a  first-order 
system  which  is  based  on  an  expanded  equation  of  motion  of  a 
structure  (3:921-922)  and  derive  the  equations  for  the  linear 
regulator  theory.  The  expanded  equation  of  motion  is  developed  from 
converting  a  structure's  finite  element  model  into  a  state-space 
system.  A  finite  element  model  which  includes  viscoelastic 
materials  described  by  a  fractional  calculus  model  will  generate  a 
state  vector  which  includes  both  integer  and  fractional  derivatives. 
The  general  form  of  a  system  created  by  expanding  a  structure's 
equation  of  motion  is  a  system  of  fractional  order  less  than  one  and 
i  s 

Dj  h  x(t)  =  A(t)  x (t)  +B(t)  u(t)  (4) 

1  Al 

where  Dt  is  a  fractional  differential  operator  of  order  1/n 
operating  on  the  state  vector  in  which  n  is  the  positive  denominator 
integer  of  the  material’s  rational  fraction  in  the  fractional 
derivative  model  of  the  material's  viscoelastic  behavior.  The  state 
vector  includes  fractional  derivatives.  It  is  shown  herein  that  this 
approach  leads  to  a  simple  control  law. 
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II.  Background 


In  the  last  ten  years  there  has  been  a  lot  of  work  on  the 
fractional  calculus  model  of  viscoelasticity.  Bagley  and  Torvik 
(1,2,3,4,15)  have  published  a  number  of  papers  on  the  model's  use 
and  accuracy.  Bagley's  PhD  dissertation  (1)  details  the  calculation 
of  material  response  for  structures  using  the  fractional  calculus 
model.  His  research  shows  that  the  model  predicts  both  oscillatory 
and  non-oscillatory  relaxation  motion  for  viscoelastic  materials 
(1:92-104).  This  latter  motion  is  due  to  the  characteristic 
relaxation  which  deformed  viscoelastic  materials  exhibit  in 
returning  to  a  less  stressed  state  over  time.  For  example,  the 
internal  structure  of  a  rubber  band  which  has  been  left  stretched  for 
a  long  time  will  slowly  adjust,  such  that  the  stress  is  relieved. 

This  has  been  observed  to  be  closely  modelled  by  power-law 
phenomenon  over  finite  observation  times  (20:681).  Furthermore, 
Bagley  (3:918)  has  experimentally  determined  the  model  parameters 
for  over  130  materials  and  has  found  the  model  accurate  over  four  or 
more  decades  in  frequency.  In  their  1985  paper  (3:921-923)  Bagley 
and  Torvik  also  illustrated  a  method  to  convert  a  finite  element 
equation  of  motion  for  a  structure  incorporating  elements  using  a 
model  of  order  a,  where  a  =  m/n  (m,  n  =  integers),  into  a  fractional 
state-space  system  of  order  1/n.  This  forms  the  basis  for  building 
the  first-order  system  required  to  develop  the  control  theory. 


Thus  far,  Skaar,  Michel,  and  Miller  (24:348)  and  Hannsgen, 
Renardy,  and  Wheeler  (12:32-34)  have  investigated  system  control 
using  the  fractional  calculus  model.  Neither  use  optimal  control 
theory  to  develop  an  optimal  control  law.  Both  assume  feedback 
forms  and  analyze  the  effects.  The  former  restrict  the  feedback 
control  law  to  a  fractional  derivative  of  the  same  order  as  the 
constitutive  law,  Pa(a)  =  P£(e)  ,  where  PCT  and  Pe  are  differential 
operators,  Pc=1  +  a3v/3tv  and  Pe=E  +  b3v/3tv ,  E  is  Young's  modulus,  a, 
b,  and  v  are  constants,  and  v  is  a  rational  fraction.  The  work  focuses 
on  the  use  of  the  root  locus  technique  to  analyze  the  system 
response.  They  derived  the  transfer  function  in  terms  of  a  new 
parameter,  h,  which  is  defined  in  the  Laplace  domain  by 
h2  =  s2PCT(s)/PE(s),  where  s  is  the  Laplace  variable.  The  poles  of  the 
transfer  function,  therefore,  are  values  of  h  instead  of  s,  the  Laplace 
variable,  and,  as  such,  cannot  be  used  directly  to  find  the  response. 
The  poles  must  be  mapped  back  to  values  of  s  through  the  definition 
of  h,  before  the  response  can  be  calculated.  A  similar  situation 
arises  in  the  application  of  the  control  theory  developed  in  this 
thesis  and  is  demonstrated  in  the  example  problems.  Furthermore, 
Skaar  also  noted  that  general  stability  criteria  can  be  defined  in  the 
complex  h-plane,  so  that  the  poles  could  be  used  directly  to  assess 
general  characteristics  of  the  system  stability.  Specificly,  Skaar 
maps  the  imaginary  axis  from  the  complex  s-plane  to  the  h-plane 
and  proves  that  all  poles  must  lie  to  the  left  of  this  mapped  line  to 
have  a  stable  system.  Similar  criteria  are  shown  in  this  thesis. 
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In  the  latter  work  referred  to  above,  Hannsgen,  Renardy,  and 
Wheeler  investigate  the  response  of  a  viscoelastic  torsional  rod 
subjected  to  linear  feedback  of  the  angular  displacement.  In  a 
special  case,  they  use  the  fractional  derivative  model  for  the 
constitutive  law  of  the  rod.  High  and  low  modes  of  vibration 
experienced  different  effects  as  control  gain  varies  from  low  to 
high  for  this  constitutive  law.  Specifically,  for  low  gains  the  low 
frequency  modes  are  least  damped;  but,  as  gain  increases,  higher 
frequency  modes  become  the  least  stable.  Optimal  gain  is  defined  as 
that  gain  which  makes  the  poles  of  the  least  stable-  mode  most 
stable  compared  to  other  gain  value  least  stable  pole  locations,  i.e. 
stability  is  classified  by  the  least  stable  mode's  most-negative-real 
location  over  the  range  of  gains  investigated.  They  define  optimum 
is  when  it  is  most  stable.  This  is  a  performance  index,  which 
optimizes  a  predetermined  control  law  instead  of  creating  one. 
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lli-  Theory 


mce  ln< 


Optimization  requires  some  measure  of  performance  to 
minimize  or  maximize.  This  performance  index  (PI)  contains 
expressions  which  represent  important  system  characteristics,  e.g. 
kinetic  energy,  and  which  are  written  in  terms  of  system  and  control 
law  parameters.  Furthermore,  these  expressions  are  weighted,  so 
that  their  effect  on  the  overall  performance  index  can  be  varied  to 
emphasize  certain  parameters  to  find  some  "best"  combination. 

To  create  a  performance  index  one  must  identify  the  important 
system  characteristics  and  their  representations  of  the  system  of 
interest.  For  this  thesis  that  system  is  a  state-space  system  of 


d| 7 n  x(t)  =  A  x ( t)  +  B  u(  t) 


where  x(t)  is  the  state  vector,  which  includes  both  integer  and 
fractional  derivatives  of  position,  and  u(  t)  is  the  control  vector.  The 
matrices  A  and  B  must  be  constant  and  real,  because  the  system 
represents  a  finite  element  model  of  a  structure.  The  structure 
equation  of  motion  has  constant,  real  coefficients.  When  the 
equation  is  expanded  into  a  state-space  system,  matrices  A  and  3 
are  formed  using  these  constant,  real  coefficients. 

Structures  are  normally  characterized  by  the  amount  of  energy 
which  the  configuration,  including  control  inputs,  stores  at  any  point 
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of  time.  Consequently,  a  performance  index  based  on  energy  would 
be  appropriate.  Current  optimal  linear  control  theory  gives  a 
starting  point.  Kirk  (15:91)  uses  a  general  PI  of  the  following  form 
for  first-order  linear  systems  of  continuous  linear  regulator 
problems: 

tf 

J  =  jxT(t|)  Hx (t[)  +  J  J (xTQx  +  UTR  u)dt 

*°  (6 
where  J  is  the  performance  index,  x(t)  is  the  state  vector,  u(t)  is 
the  control  vector,  and  H,  Q,  and  R  are  real,  symmetric  weighting 
matrices,  in  which  R  is  positive  definite  and  H  and  O  are  positive 
semi-definite.  Values  of  H,  Q,  and  R  are  chosen  for  an  optimization: 
consequently,  each  set  of  weighting  matrices  defines  an  optimized 
value  of  the  performance  index,  J.  For  this  thesis  H,  Q,  and  R  are 
chosen  to  be  constant  matrices.  Since  the  state  vector  contains 
both  position  and  its  first  derivative,  1/2  xT  Qx  includes  measures 
for  both  the  strain  and  kinetic  energies  of  the  structure.  Similarly, 
1/2  uT  Ru  is  a  measure  of  energy  in  the  controller.  The  expression 
with  H  merely  represents  the  energy  in  the  desired  final  state. 
Although  the  extra  fractional  states  are  now  included  in  the  system, 
the  classical  definitions  of  strain  and  kinetic  energy  are  assumed  to 
still  be  valid.  Furthermore,  the  extra  fractional  states  are  not 
known  to  represent  any  energy  storage  system.  That  is,  the 
structure  stored  energy  is  still  described  by  kinetic  and  strain 
energies.  Consequently,  a  quadratic  performance  index  should  still 
be  suitable  for  optimizing  the  1/nth-order  system.  Weighting 


1  0 


matrices  H  and  O,  however,  must  be  chosen  such  that  the  fractional 
states  have  zero  weight.  This  eliminates  these  states  from  the 
index. 

For  the  problems  of  interest,  a  slightly  modified  form  of  this 
PI  is  desireable.  For  simplicity  the  transient  solution  of  the  state 

vector  is  not  addressed,  and  only  the  steady-state  solution  is 
sought.  Therefore,  let  time  be  large  and  let  tf  tend  to  infinity. 

Furthermore,  when  the  desired  final  state  is  static  equilibrium  of 
the  structure,  1/2  x"*"(tf)  Hx(tf)  is  assumed  to  be  equal  to  zero. 
Therefore,  the  actual  performance  index  used  herein  to  optimize 
fractional-order  systems  of  the  form  of  Eq  (5)  is 


oo 


(7) 


Kirk  (15:86-90)  thoroughly  discusses  the  optimization  of  his 
general  PI,  which  is  a  more  general  form  of  the  chosen  performance 
index, 


J(x(t),  t,  u(t))  =  h(x(tf),  tf)  +  g(x(x),  u(t),t)  dt 

t<TStf  Ji 


for  a  general  first-order  system  of  the  form 


x(t)  =  a(x(t),u(t),  t) 

If  the  1 /nth-order  system  can  be  posed  as  a  first-order  system,  then 
the  optimization  process  for  our  performance  index  has  already  been 
outlined  by  Kirk.  The  development  can  be  followed,  and  equations 


pertaining  to  the  system  of  order  1/n  can  be  found.  The  linear 
regulator  problem  is  specificly  followed  because  it  uses  the  same 
performance  index  as  has  been  chosen  for  this  thesis  and  results  in 
linear  feeback  of  the  state  matrix.  It  will  be  shown  that  the 
modified  control  law  will  also  be  linear  or  asymptotically  linear  for 
steady-state  (t  -*«>). 

Once  a  structure  has  been  built,  the  only  way  to  affect  its 
stability  is  through  applying  some  control  force  or  torque,  u(t).  The 
state  vector,  x(t),  is  just  a  statement  of  the  location  and  motion  of 
the  structure's  parts,  i.e.  it  is  a  description  of  the  state  of  the 
structure.  The  control  u(t)  is  what  is  being  implemented  to  cause  a 
change  in  the  state  of  the  structure.  It  is  also  the  quantity  for 
which  the  physical  control  system  is  going  to  be  designed. 

Therefore,  the  control  law  which  minimizes  the  performance  index 
is  what  is  being  sought.  Consequently,  the  performance  index  is 
minimized  with  respect  to  the  control,  u(t).  The  PI  is  minimized, 
because  the  control  objective  is  minimum  structural  motion  and 
minimum  energy  used  by  the  controllers,  all  of  which  is 
characterized  by  minimum  energy  in  the  system.  The  control  u(t) 
which  minimizes  the  index  is  designated  the  optimal  control,  u*(t). 
The  relation  between  u*(t)  and  its  variables  (e.g.  x(t)  and  t)  is  called 
the  optimal  control  law.  Therefore,  the  basic  concept  behind  the 
development  of  this  optimal  control  theory  for  systems  of  order  1/n 
is  the  determination  of  an  optimal  control  law  by  minimization  of 
the  performance  index  with  respect  to  the  control,  u(t): 


1  2 


(9) 


J  (x(t),t) 


aw  « + 

t<T<tf  ' 


g(x(x),  u(x),x)  dx 


rtf 

-  h(x(Df),  tf)  +J  g(x(x),  u  (x),  x)  dx 

* 

where  J  (x(t),t)  is  the  optima!  performance  index  function.  It  is  no 
longer  a  function  of  u(t),  since  the  minimum  over  all  admissible 
controls  has  been  taken;  i.e.  J*  specifies  u*(t)  as  the  optimal 
control. 

*  _* 

Eq  (9)  defines  J  (x(t),t).  However,  the  equation  includes  two 

unknown  functions,  because  the  optimal  performance  index  function 

is  unknown  as  well  as  the  optimal  control.  There  is,  thus  far,  one 

equation  for  two  unknowns;  consequently,  an  expression  for  one  of 

*  _ # 

the  unknowns  must  be  assumed.  A  function  is  chosen  for  J  (x(t),t)  to 
emphasize  system  characteristics  whcih  are  believed  to  be 
important,  and  the  corresponding  u*(t)  is  found. 

*  m 

If,  instead,  the  minimization  is  performed  with  J  (x(t),t) 

*  t  r 

unspecified,  the  optimal  control  will  be  a  function  of  J  (x(t),t),  x  ( t), 

* 

and  t.  Depending  on  its  form,  J  (x(t),t)  may  have  one  or  more 
constants  of  proportionality  between  itself  and  its  variables; 
therefore,  one  or  more  equations  are  still  necessary  to  completely 
describe  the  optimal  control.  If  a  simple  form  for  the  optimal  index 
function  is  assumed,  e.g.  J  (x(t),t)  =  1/2  xTR(t)x  ,  then  one  equation 
for  K(t)  and  one  for  u*(t),  including  K(t)  and  x(t),  are  required.  This 
is  the  case  for  th9  linear  regulator  problem.  K(t)  is  often  called  the 
optimal  "gain"  matrix,  because  it  is  the  matrix  of  gains,  i.e. 


amplification  factors  applied  to  a  feedback  state,  for  the  assumed 
optimal  performance  index  function.  Both  equations  are  found 
through  the  minimization  of  the  index  as  shown  in  the  next  section. 

Kirk  illustrates  the  optimization  of  a  general  performance 
index  with  the  general  first-order  system,  Eq  (9).  If  there  is  an 
equivalent  first-order  system  for  the  1  /nth-order  system,  it  can  be 
used  the  general  equations  that  result  from  the  optimization.  Since 
the  optimization  is  simple  and  Kirk  (15:86-87)  describes  it  well,  it 
is  not  repeated  here.  Instead,  this  thesis  will  start  from  the 
resulting  equation,  the  Hamilton-Jacobi-Bellman  equation,  and  will 
detail  the  modification  of  the  linear  regulator  theory  equations  for  a 
system  of  order  1/n. 

Modified  Hamilton-Jacobi-Bellman  Equation  and  Optimal  Control  Law 
Kirk's  development  leads  to  the  important  Hamilton-Jacobi- 
Bellman  (H-J-B)  equation,  which  must  be  satisfied  by  the  optimal  PI, 
J*(x(t),  t): 

0  =  jf  (x(t),  t)  +H(x(t),  u  (x(t),  Jx,  t),  Jx,  t)  (10) 

where  the  assumed  final  value  is  J*(x($),  ^)  «  0  and  H ,  called  the 
Hamiltonian,  is 

H  (  x(t),  u(t),  Jx  ,  t)  =  g(  x  ( t) ,  u(t),  t)  +  Jx  (  x  ( t) ,  t)[a  (  x  ( t) ,  u ( t) ,  t )]  (i  i) 

where  a(x(t),u(t),  t )  is  a  first-order  system.  The  Hamiltonian  is 
minimized  by  the  optimal  control,  u*(t): 


Eq(10)  and  Eq  (11a)  are  important,  because  the  minimization  of  the 
Hamiltonian  yields  the  optimal  control  law  and  the  H-J-B  equation 
yields  an  equation  for  the  optimal  gain  matrix  in  the  performance 
index.  Hence,  these  are  the  two  equations  sought  in  the  previous 
section  to  describe  the  optimal  control  law. 

The  H-J-B  equation  is  the  relation  that  defines  the  general 
optimal  performance  index,  J*(x(t),  t).  However,  as  noted  earlier,  a 
form  must  be  assumed  for  J*(x(t),  t),  e.g.  1/2  xT  K(t)  x  where  K(t)  is 
the  optimal  gain  matrix.  If  the  PI  is  chosen  appropriately,  the  H-J-B 
equation  will  give  an  equation  for  K(t).  Kirk  pursues  this  approach. 
First,  he  minimizes  the  Hamiltonian  with  respect  to  a  general  J*  and 
finds  that  the  optimal  control  law  which  minimizes  it  is  given  by 
linear  feedback  of  the  partial  derivative  of  J*  with  respect  to  x(t). 
Thus,  the  linear  regulator  problem  is  defined.  If  J*(x(t),  t)  is  chosen 
to  be  equal  to  1/2  xTK(t)  x,  the  optimal  control  is  linear  feedback 
of  the  state  vector.  Using  this  optimal  control  and  PI,  the  H-J-B 
equation  leads  to  a  nonlinear  ordinary  differential  equation  for  the 
optimal  gain  matrix,  K(t).  For  the  linear  regulator  problem  this 
equation  is  a  Riccati  equation.  If  the  1/nth-order  system  could  be 
posed  as  a  first-order  system,  then  a  development  parallel  to  Kirk's 
can  be  used  to  define  modified  versions  of  the  H-J-B  equation, 
optimal  control  law,  and  Riccati  equation.  Therefore,  the  problem 
reduces  to  writing  the  1  /nth-order  system  in  the  form  of  a  first- 


order  system.  Then  substitute  this  expression  and  the  index  form, 

Eq  (6),  chosen  for  this  thesis,  into  the  general  Riccati  and 
Hamiltonian  equations,  and  complete  the  minimization  for  this 
specific  linear  regulator  problem. 

The  first-order  system  can  simply  be  rewritten  recursively  in 
terms  of  the  derivative  of  order  1/n  which  was  defined  in  Eq  (1): 

x  =  d{  *  [d{  *  ...  (dJ  *(x)  )]  (n  times)  (12) 


Successive  replacement  of  D{  ^x  by  A  x(t)  +Bu(t),  as  in  Eq  (4), 
and  application  of  the  next  derivative  operator  yields  a  general 
first-order  system  of  the  form: 


where  D°  u(t)  =  u(t).  This  is  illustrated  by  an  example  of  a 
1 /3-order  system  with  constant  matrices  A  and  B: 

-  -.1X3  f-,1  X3  /—.I  /3/  — \  ) 

x  -  Dt  [Dt  \Dt  (xj/j 
x  =  D1,  a  [d{  a  (  A  x{t)  +B  u(t)  ). 
x  =  D1,  73  [  A  D1,  13  x  (t)  +  B  D{  13  u(  t)  ] 
x  =  D{  a  [  A  (A  x  ( t)  +  B  u(  t))  +  B  D{  a  u(  t) 


(1 4.a-g) 


x  =  A2  d{  a  x  (t)  +  AB  d{  a  u(  t)  +  B  Dt  13  u(  t) 
x  =  A2  (  A  x(t)  +  B  u(  t)  )  +  AB  D{  a  u(  t)  +BC?/3u(t) 
x  =  A3  x  (t)  +  A2  B  u(  t)  +  AB  D{  **  u(  1)  +BD^/3  u(  t) 


The  resulting  expression  can  be  used  with  the  H-J-B  equation  to 
develop  the  optimal  control  law  and  a  modified  Riccati  equation 
governing  the  gain  K(t)  in  the  1/nth  order  system.  An  important 
point  is  that  this  development  is  restricted  to  systems  of  rational 
1 /nth-order  since  n  must  be  an  integer  to  apply  the  derivative  n 
times  in  the  recursive  relation  for  x. 

The  first  step  in  using  this  equivalent  form  in  the  H-J-B 
equation  is  minimization  of  the  Hamiltonian  with  respect  to  the 

control  vector;  a  necessary  condition  is  vanishing  of  the  first 

dH  d2H 

variation:  rr  =  0  .  If  the  second  variation,  r  ,  is  a  positive- 

3j2 

definite  quadratic  form,  the  control  vector  which  solves  this  is  the 
desired  optimal  control  vector.  The  minimum  value  of  H  is  found  by 
evaluating  it  at  the  optimal  control  vector. 

Using  the  equivalent  first-order  system,  Eq  (13),  the 
Hamiltonian  becomes 


H  =  J-xtQx 


-  *T 

Ru  +Jx 


n-  1 

Ax  +  j 
p=o 


-n-p  -  1 


b  dj 


/  ri¬ 
ll 


(15) 


To  find  the  first  variation  with  respect  to  ij  ,  evaluate  H  at  u  +  h  : 


1  7 


H  =  l-xTGk  +l(u  +h)TR(u  +  h)  +  Jx  (  Anx  +  J  Anp  1 BD^7 n(u  +h 
22  \  p=o 


— T' 


Simplify,  using  1-  h  Ru  =^-uTRh  (which  is  a  scalar);  the  result  is 
2  2 


H  =  H  (u)  +  —  (h)  +  —(h) 


8u  9u2 


where 


(16a) 


H  (u)  =  J-xtO<  +—  utRu  -kIx  [  Anx  +  2^  An  P  1  Bdf  nu 
22  V  p=o 


dH 

dU 


^  *T  n- 1 


(h)  =  u TRh  +  Jx  j  An'p'1BC^/nh 


p=0 


(16b) 


(16c) 


™(h)-ihTK 

au2  2 


(1 6d) 


This  is  an  the  expansion  of  H  evaluated  at  u  +  h.  The  first  variation, 
or  derivative,  of  H  operating  on  h  corresponds  to  the  terms  linear  in 
h  .  The  second  variation  corresponds  to  the  terms  quadratic  in  h  . 
With  the  requirement  that  R  be  positive  definite,  the  second 
variation  is  always  positive.  Therefore,  a  solution  of  the  first 
variation  will  give  a  minimum  value. 

The  first  variation,  Eq  (16c),  operating  on  h  must  equal  zero 
for  any  admissible  h  , 
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where  admissible  controls  and  variations  are  taken  to  be  functions 
h(t)  which  are  continuous  for  0  <t<°«.  Since  Eq(17)  must  be  true 

for  all  admissible  variations,  h(t)  ,  it  must,  in  particular,  be  true  for 
the  admissible  variation 

h(t)  =  H(t)  T  (18) 


where  H(t)  is  the  Heaviside  step  function  which  is  zero  for  t  <  0  and 
1  for  t  >  0;  the  vector  1  is  the  constant  vector  having  the  value  1  for 
all  components.  For  this  admissible  variation,  the  fractional 
derivative  is 


□P'n  h<  t)  .  rf,nH(l)  7- - 3 - jd_  (  H(  frt)  t'P' n  dt  1 

r(  1  -p/  n)  dt  J0 

which  yields 

D?'nh(t)-Tp„(t)  7 
where 

I  0 ,  ISO 

Tpn  ( t)  =  )  1,  P=0,t>0 

\ — —  ,  p  >  0  ,  t>  0 

Ir(1-p/n) 


(19) 


(20) 


Using  Eq  (20)  in  Eq  (17),  the  optimal  control  u  must  be 


(21) 


(22) 


u*(t)  =  -  R  '  I  BT(An'P'')TTpn(t)  i’(  x(t),  t, 


p=0 


Observe  that  the  large-time  behavior  of  u  is 


u*(t)~  -R  1  BT(a  n  ^  Jx(  x(t),  t)  ,  t-> 


(23) 


The  small-time  form  of  u  includes  singular  factors  through  Tpn(t) 

-  * 

which  must  be  dominated  by  appropriate  behavior  of  Jx  for  small 
values  of  t,  in  order  for  u  to  be  an  admissible  control.  This 
constraint  on  J  will  appear  in  initial  conditions  on  the  gains  K 
obtained  from  a  Riccati  equation. 


IM  Riccati  Equation  for  the.  Gains 


The  goal  of  this  section  is  to  derive  the  Riccati  equation  which 
governs,  for  large  time,  the  gain  matrix,  K(t),  for  linear  regulator 
problems  in  which  gains  are  asymptotically  constant.  Upon  choosing 
an  appropriate  form  for  the  optimal  PI  and  applying  it  to  the  control 
law,  Eq  (22),  the  Riccati  equation  is  derived  with  an  application  of 
the  H-J-B  equation.  Using  Eq(16b)  in  Eq  (9),  the  H-J-B  equation  is 
now 


0  =  4 


x  TCfc  +  —  u 
2  2 


*T * 

Ru 


+Jx 


~n  H-  1 

A  Z  +  t 

p=0 


'n-p  - 1 


b  dj 


/  n  — 
U 


\ 

(24) 


Recall,  linear  state  feedback  is  achieved  if  J*  is  quadratic  in  x. 
Furthermore,  it  is  logical  that  the  PI  may  be  related  to  the 
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expressions  representing  the  system  characteristics  of  interest. 
Therefore,  it  is  commonly  assumed  that  the  PI  takes  the  form 


J  ( x  ( t) ,  t)  =  x T  K  ( t)  x 


(25) 


where  K(t)  is  real,  symmetric,  and  positive-definite.  Then 


Jx  =  Kx 


(26) 


jJ  =  J-xTKx 
2 


(27) 


Using  Eq  (26)  in  Eq  (22),  the  optimal  control  law  becomes 


Pi  *1  r->Tf  a  n'P*  1IT 


u  ( t)  »  -  R  '  X  B  (A 
p=o 


Tpn(t)  K(t)  X  (t) 


(28) 


(From  this  result,  the  large-time  behavior  of  u  can  be  viewed  in 
terms  of  the  behavior  of  x,  since  K  is  asymptotically  constant. 
Furthermore,  the  initial  behavior  of  K  can  be  ascertained  to 
dominate  the  singular  nature  of  Tpn(t).) 

When  Eq  (27)  and  the  optimal  control,  u  ,  in  Eq  (28)  are 
inserted  into  the  H-J-B  equation,  Eq  (24),  the  result  after 
simplification  is 

ixT{i<-fQ+R2  An'p''  BTpn(t)  R  1 BT  X  (a"'p'1  )TT,„(t)  K  +  2KA" )  x 

d  \  p-o  q-0  / 
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n-1 


^xT  2K  X  An'p*1B  Df/n 

2  p-0 


n-1 


I  R'1BT(Ar"q‘1)  Tqn(t)  K(t)x(t) 
L  q-o 


(29) 


At  this  point,  we  note  that  a  commonly-employed  technique  saying 
"if  x  tMx  =  0  for  any  x,  then  M  =  0"  cannot  be  used  since  the 
differential  operator  D^/n  operates  on  both  xand  K;  they  cannot  be 
separated  in  Eq  (29)  to  use  this  argument  to  get  an  equation  which  K 
must  satisfy.  However,  consideration  of  the  individual  terms 
summed  over  p  in  the  right  side  of  Eq  (29)  reveals  that,  for  gains  K 
which  are  asymptotically  constant  for  large  time,  there  are  states, 
x(t),  for  which  the  leading  order  behavior  of  this  sum  is  given  by  its 
p=0  term.  In  effect,  the  fractional  derivatives  of  order  greater  than 
zero  are  subdominant  to  the  undifferentiated  (p=0)  summation  over 
q.  Relying  on  "regularly-varying"  function  theory  (Appendix  A), 
Warhola  (27)  illustrates  states  for  which  this  is  true,  so  that  the 
theory  presented  hereafter  is  not  empty.  Proceeding  with  this 
Eq  (29)  yields 
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\  xT  j  K  +  Q  +  K  X  An'p'1  BTpn(t)  R’1  §T  £  (a"'"''  )TT,„(t)  K  +  2KAn )  x 

d  (  P-0  q-o  ) 


^  lxT2KAn‘1B  £  R‘1BT(An'q‘1)TTqn(t)Kx  ,  t  -»  oo 
*  q=o 


(30) 


The  leading  order  of  the  sum  over  q  is  clearly  the  zeroth  term,  so 
4-xT  j  K  +  Q  +  K  £  An’p'1BTp„{t)  R-'bT  £  (An'<,'1)TTqn(t)  K  +  2Ra")  x 

d  (  P-0  q-0  ) 


xT  2K  An1  B  R  1 BT  ( A0*1  )T  K  x  ,  t -»  00 


(31) 


Now  this  must  be  true  for  all  states,  x(t),  so 


K  +  Q  +  K  X  An'p'1BTpn(t)  R 1BT  £  (An‘q'1)TTqn(t)  K  +  2KAn 

p»0  q-0 


^  2KAn_1B  R1BT(An'1)TK  ,t 


— >  00 


(32) 


Similar  examination  of  the  terms  in  the  left  side  of  Eq  (32)  reveals 
that  the  p=q=0  term  dominates  for  large  time,  at  least  for 
asymptotically  constant  gains,  K,  which  have  a  regularly-varying 

«N> 

derivative,  K.  (In  fact,  the  K  term  is  subdominant  to  all  terms  with 
p+q  <  n  on  the  left  side  of  Eq  (32).)  Thus, 
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Q+K  A^B  R‘1BT(An'1)TK  +  2KAn 


'n-1~  ~ -1~T  f ~n- 1 


^  2K  A  B  R  B  (A  K  ,  t->  oo 


(33) 


In  the  limit  as  t  ->  oo  ,  the  algebraic  Ricca 
steadv-state  gains.  K,  is  obtained  as 

X  ~Tn  ~  rn-1~  ~-1~T/rn-l\T~ 

Q+  2KA  -KA  BR  B  |A  )  K  -  0 


ling  the 


(34) 


Now,  the  matrix  KA  can  be  split  into  its  symmetric  and  anti¬ 
symmetric  parts 

KAn  _  J-  [kA n  +(KAn)T]+l[KAn  -  (ka"}1] 

2  2 

Using  the  fact  that  the  transpose  of  a  scalar  is  itself,  i.e. 
xtk£"x  =  xT(KAn)Tx 


then  2KA  can  be  written  as 


2KAn  =  KAn  +  (KAn)T 

Substituting  Eq  (35)  into  Eq  (34)  and,  again,  using  the  fact  that  the 
transpose  of  a  scalar  is  itself,  the  result  is  the  final  form  for  the 


>raic  Hiccati 


tion  oovermr 


0=  Q+KAn +(KAn)T  -KAn'1BR'1B1|An‘1)TK 


(36) 
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It  is  important  to  note  that  the  constant-K  Riccati  equation 


for  first-order  systems  is 

0  =  K  +  Q+KAi  +(KA1f  -KBtR^bIk  (37) 

•  where 

x(t)  =  Ai  x(t)  +Bi  u(t)  (38) 

Consequently,  the  modified  1/nth-order  system  equation  can  be 

•  written  as  ..an  equivalent  first-order  system  equation,  if 

Ai  =  An  (39) 

B1-An'1B.  (40) 


The  significance  of  this  result  is  that  current  algebraic 
Riccati  equation  solving  routines  for  first-order  systems  can  be 
used  to  solve  the  modified,  constant-K,  steady-state  equation  for 
1 /nth-order  systems.  The  tools  to  do  this  analysis  are  already 
available  and  in  common  use. 


25 


Solution  Process  for  System  Response 


Once  K  has  been  found  by  solving  the  Riccati  equation  (36),  the 
solution  for  the  state  vector,  x(t),  can  be  completed.  This  solution 
is  typically  called  the  response  of  a  system,  because  the  state 
vector  describes  how  a  structure  reponds  to  a  given  input.  By 
substitution  of  the  steady-state  optimal  control  expression, 

Eq  (29),  into  the  system  of  order  1/n,  Eq  (5),  the  steady-state 
system  equation  becomes 


Dt  x  = 


A  -  R  BA 


K 


(41a) 


The  Laplace  transform  of  Eq  (41a),  assuming  zero  initial  conditions, 
i  s 


»i  /ft 


X(s)  = 


A  -R^bIa"'1)  K 


X(s) 


(41b) 


where  s  is  the  Laplace  variable  and  capital  letters  on  the  vectors 
indicate  transforms.  In  the  structure  examples,  the  state  vector  is 
actually  a  composite  of  fractional  order  derivatives  of  displacement 
vectors,  e.g. 


where  (x)  is  a  vector  of  displacements.  The  solution  technique  is 
discussed  next  in  general  terms  and  is  illustrated  in  detail  in  the 
example  problems  in  Chapter  4. 

If  the  system  eigenvalues,  X,  are  defined  by  X  =  s1/n  ,  then  the 
solution  process  is  reduced  to  an  eigenvalue  problem,  i.e.  Xx  =  Cx  , 
where  C  is  a  real  matrix  for  structures: 


XX(s)  = 


r  s-isurn-i 

A  -  R  BA 


Kj  X(s) 


(42) 


Insight  to  solving  Eq  (42)  for  the  state  vector  response  can  be 
gained  by  reviewing  the  eigenvalue  problem  of  a  first-order  system. 

In  the  classical  first-order  system  eigenvalue  problem 
eigenvalues  are  defined  in  the  Laplace  domain  by  the  relation  co  =  s  , 
which  results  in  the  solution  of  the  state  vector  as  x  =  e cot  , 
where  <t>  is  a  vector  of  coefficents  and  co  is  a  complex  number.  The 
vector  <t>  is  called  an  eigenvector,  and  each  eigenvalue  defines  one. 
The  complex  plane  in  which  the  eigenvalues  lie  will  be  called  the 
s1 -plane,  where  the  superscript  indicates  the  power  of  the  Laplace 
variable  s  in  the  definition  of  the  eigenvalue.  If  Eq  (42)  were  an 
eigenvalue  problem  of  a  first-order  system,  by  substitution  of  the 
transform  of  the  general  solution  of  the  state  vector,  Eq  (42)  can  be 
written  as 


0  = 


~  -1~Tf~n- 1  'T~ 


A  -  R  B  \A  )  K  -  X\\<t> 


(43) 
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« 


where  I  is  the  identity  matrix.  If  the  eigenvector  is  not  the  trivial 
solution,  Eq  (43)  has  a  solution  if  and  only  if  the  determinant  of  the 
matrix  of  coefficents  of  the  eigenvector  is  zero  (18:78). 


0  = 


~  ~-1~Trn-1 

A  -  R  BA 


K  -  XI 


(44) 


Eq  (44)  yields  a  polynomial  equation  in  the  eigenvalue  X,  and  the 
roots  of  this  polynomial  are  the  eigenvalues.  Because  all 
coefficients  of  the  polynomial  are  real,  the  eigenvalues  may  include 
real  values  or  pairs  of  complex  conjugate  numbers.  Once  the 
eigenvalues  are  determined,  the  eigenvectors,  <!>,  can  be  determined 
from  Eq  (43).  The  vectors  are  unique  in  that  the  ratio  between  any 
two  elements  is  constant,  but  the  values  are  not  unique  (18:78).  The 
ratio  is  dependent  on  the  eigenvalue,  but  the  values  are  arbitrary 
because  c  <t>,  where  c  is  an  arbitrary  constant,  is  also  a  solution  of 
Eq  (43).  In  typical  equations  of  motion  for  structures  x,  hence  <t>,  are 
displacements.  Since  only  the  ratio  of  elements  in  the  eigenvector 
is  constant,  the  eigenvector  represents  relative  displacements 
between  nodal  ponts  in  the  structure.  If  the  eigenvector  is  real 
valued,  this  gives  a  shape  of  the  structure.  Consequently,  the 
eigenvector  is  called  a  mode  shape,  i.e.  the  shape  of  a  mode  of 
motion.  These  modes  are  quantized  in  that  they  only  occur  in 
conjunction  with  specific  eigenvalues.  It  can  be  shown  that  the 
eigenvalues  correspond  to  frequencies  of  vibration,  if  the  motion  is 
oscillatory  (complex  eigenvalue),  or  reciprocals  of  time  constants, 
if  the  motion  is  non-oscillatory  (real  eigenvalues)  (18:78). 
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Therefore,  structures  may  vibrate  in  a  specific  shape  at  one  of  these 
"natural"  frequencies  or  have  a  specific  shape  as  the  motion  decays 
or  grows  with  a  specific  time  constant.  Modes  are  named  by  the 
type  of  motion  and  eigenvalue  ranking,  when  eigenvalues  are  ordered 
starting  with  the  lowest  value  for  that  type,  e.g.  the  second 
vibrational  mode  corresponds  to  the  next  eigenvalue  higher  than  the 
lowest  vibrational  eigenvalue.  If  the  eigenvector  is  complex,  the 
mode  shape  changes  as  time  passes.  Finally,  the  total  structural 
response  may  be  expressed  as  a  linear  combination  of  responses  due 
to  individual  modes  (19:164).  The  lowest  modes  have  the  greatest 
influence  because  displacements,  hence  energy,  are  greater  than  the 
higher  numbered  modes.  Once  eigenvalues  and  eigenvectors  are 
known  in  the  classical  eigenvalue  problem,  one  can  go  directly  to 
response. 

Consequently,  eigenvalues,  called  poles  in  control  analysis,  and 
eigenvectors  are  typically  the  focus  in  structural  and  control 
analysis.  Because  the  pole  affects  the  exponent  of  the  exponential 
in  the  solution,  x  =  pe®1  .its  value  determines  stability,  i.e. 
bounded  motion,  of  the  response.  If  co  is  complex,  the  exponential 
can  be  written  as  a  product  of  two  exponentials:  one  with  a  real 
exponent  and  one  with  an  imaginary  exponent.  The  Euler  formula 
(17:584)  states  that  the  latter  can  be  written  in  terms  of  a  real 
cosine  part  and  an  imaginary  sine  part.  Therefore,  a  complex  pole 
produces  oscillatory  response,  in  which  the  magnitude  of  the 
imaginary  coefficient  of  time,  t,  in  the  exponent  is  the  frequency  of 
oscillation.  Oscillatory  motion  is  bounded,  so  it  is  the  real 
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component  of  the  pole  which  must  determine  the  boundedness  of 
motion.  If  the  real  part  is  positive,  its  exponential,  which  is  the 
amplitude  of  the  oscillatory  motion,  grows  without  bound.  If  the 
real  part  is  negative  or  zero,  the  response  is  bounded,  and  the 
system  is  stable.  Therefore,  pole  location  can  give  a  general  idea  of 
the  response  of  the  function,  e.g.  a  pole  in  the  second  quadrant  is 
stable  and  produces  oscillatory  motion,  and  is  a  primary  tool  in 
analysis  of  structural  motion  and  control.  The  criteria  that  non¬ 
positive  real  parts  of  poles  produce  bounded  motion  only  applies 
when  the  poles  are  expressed  in  this  plane.  Other  criteria  will  be 
derived  for  eigenvalues  expressed  in  other  planes. 

For  systems  of  order  1/n,  the  eigenvalue  is  defined  herein  by 
X  =  s1/n  ,  where  X  is  the  eigenvalue  and  the  eigenvalues  lie  on  a 
complex  plane  called  the  s1/n-plane.  The  general  form  of  the  time- 
domain  state  vector  response  is  no  longer  the  exponential;  it  is 
something  more  complex.  A  solution  to  Eq  (42)  can  still  be  found 
relatively  easily  by  mapping  the  eigenvalues  for  the  system  of  order 
1/n  to  a  first-order  system.  Recalling  the  recursive  equation  for  an 
equivalent  first-order  system,  Eq  (12),  one  can  determine  that  for 
an  eigenvalue,  X,  which  is  defined  by  first-order  system,  X  =  con  . 
This  implies  two  important  ideas.  First,  one  can  determine  the 
eigenvalues,  X,  map  them  to  co,  and  determine  system  response  from 
co  as  described  for  the  classical  system.  Second,  because  the 
eigenvalue  X  is  the  nth-root  of  the  eigenvalue  co,  there  are  n 
^.-eigenvalues  corresponding  to  a  single  co-eigenvalue.  These  X-poles 
actually  produce  more  than  just  the  response  described  by  the 
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co-eigenvalue.  Some  of  the  A.-poles  produce  the  exponential  response 
desribed  by  a  corresponding  co-eigenvalue,  but  the  others  have  been 
shown  to  generate  terms  of  negative  powers  of  time  in  the  response 
(1:92-104).  This  power  response  correspunds  to  the  relaxation  or 
creep  response  observed  by  Nutting  in  viscoelastic  materials  (2:746, 
20:681).  If  the  equations  of  motion  for  the  system  of  order  1/n  are 
transformed  into  the  Laplace  domain,  the  complex  plane  will  be 
called  herein  the  s1/n-plane,  where  the  superscript  denotes  the 
order  of  the  derivative  in  the  eigenvalue  definition. 

The  type  of  response  due  to  a  pole  in  the  s1/n-plane  can  be 
determined  from  pole  location,  just  as  in  the  s1 -plane.  Because  the 
X-poles  will  be  mapped  into  the  s1 -plane  to  determine  response,  it 
is  important  to  recognize  which  regions  of  the  s1/n-plane  map  to  the 
four  quadrants  of  the  s1 -plane  (Figure  1). 


s^plane  in  s1/n-plane 


Figure  1. 


Mapping  Regions  Between  s1/n-plane  and  s^plane 


s 


1  -plane 
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The  type  of  response  due  to  poles  in  these  regions  of  the  s1/n-plane 
is  determined  by  the  type  in  the  corresponding  s1 -plane  quadrant. 
Poles  will  also  lie  outside  these  regions,  and  it  has  been  shown  that 
these  poles  produce  relaxation  response  described  by  the  terms  of 
negative  powers  of  time  (1:92-104).  The  next  section  addresses 
poles  and  mapping  in  more  detail. 

Mapping  Regions 

Starting  from  the  s1/n-plane  and  mapping  to  the  s1 -plane,  the 
mapping  function  is 

zi  =  (zi/n)n  where  zi/n  =  complex  number  in  s1/n-plane 

zi  =  complex  number  in  s1  -plane 

If  one  considers  the  inverse  mapping,  i.e.  z-|/n=(zi)1/n,  it  is  evident 
that  each  quadrant  of  the  s1 -plane  maps  to  a  region  of  less  angular 
width  in  the  s1/n-plane.  The  kth  quadrant,  (k-1)jt/2  <  arg(zi)  <  kjc/2 
for  k=1  to  4,  maps  to  the  smaller  sector  (k-1)7i/2n  <arg(zi)  < 
k7i/2n.  It  is  easy  to  identify  which  sectors  of  the  s1/n-plane  map  to 
each  of  the  quadrants  in  the  s1 -plane.  The  rays  in  the  s1/n-plane  and 
the  s1 -plane  axes,  to  which  they  correspond  (or  map),  bound  the 
sectors  in  their  respective  plane: 
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s^n-plane 

maps 

s^plane 

The  ray  on  which  s1/n  =  A1/n  e  i0 

-» 

to 

maps 

+  real  axis 

The  ray  on  which  s1/n  =  A1/n  e  ±  i7C/n 

— > 

to 

maps 

-  real  axis 

The  ray  on  which  s1/n  =  A1/n  e  iTC/2n 

— » 

to 

maps 

+  imaginary  axis 

The  ray  on  which  s1/n  =  A1/n  e'i7C/2n 

— » 

to 

-  imaginary  axis 

where  A  is  an  arbitrary  magnitude.  This  mapping  is  shown 
graphically  in  Figure  1.  It  is  important  to  recognize  that  regions  do 
not  reduce  or  increase  in  size.  The  regions  just  appear  differently 
on  the  two  planes.  The  s1/n-plane  view  only  shows  more  regions  so 
that  the  relaxation  response  can  also  be  seen. 

To  understand  to  what  these  other  regions  refer,  it  is 
necessary  to  discuss  the  concept  of  a  Riemann  surface  and  sheets  of 
complex  numbers  (17:621).  Complex  numbers  be  considered  as  on  a 
spiral  surface  like  an  infinite  corkscrew,  because  the  angular 
parameter  of  a  complex  number  is  not  unique.  The  angle  is  the  same 
if  2rtp  radians,  where  p  is  an  integer,  are  added  or  subtracted.  Each 
complete  twist  of  the  screw  represents  one  cycle  of  2ti  radians. 

This  twist  is  a  sheet,  often  called  a  Riemann  sheet,  of  complex 
numbers.  Where  the  sheet  begins  and  ends  depends  on  the  branch  cut, 
which  is  defined  by  the  specific  problem  and  function  singularities, 
in  the  plane.  For  example,  since  in  this  thesis  the  branch  cut  is 
assumed  along  the  negative  real  axis,  a  sheet  may  be  defined  by 
those  complex  numbers  whose  argument  (angle)  lies  between  k  and 


.K.  An  analytic  complex  function  requires  a  closed  Riemann  surface 
in  order  to  satisfy  the  requirement  of  differentiability  everywhere. 
The  surface  is  composed  of  sheets  joined  together  at  the  branch 
cuts,  such  that  a  continuous  surface  is  formed.  Because  the  function 
w  =  z  n  maps  one  sheet  to  n  sheets,  i.e.  numbers  with  argument  2n 
become  numbers  with  argument  2nn,  the  corresponding  Riemann 
surface  consists  of  n  sheets  linked  together  at  the  branch  cuts  like 
a  chain.  The  first  and  last  cuts  join  to  make  a  continuous  surface, 
like  a  loop  of  chain. 

In  the  classical  first-order  system  eigenvalue  problem  in 
which  the  eigenvalue,  X,  is  defined  by  d]x  =  Xx  and  there  are  no 
fractional  powers  of  X,  the  mapping  function  is  w  =  z1 ,  and  there  is 
only  one  sheet.  Hence,  the  reader  may  not  be  familiar  with  poles  on 
other  Riemann  sheets.  Once  the  eigenvalue  problem  is  defined  by 
d{  ^  x  =  Xx  ,  there  are  n  sheets  on  whichthere  are  have  poles;  the 
s1/n-plane  view  shows  them  all  as  sectors  of  angular  width  2rc/n 
radians.  Considering  the  sheets  as  a  stack,  it  is  easy  to  visualize 
the  plane  representation  as  the  projection  of  all  the  poles  in  the 
stack  onto  a  planar  surface.  Often,  because  one  is  interested  in 
poles  on  only  one  sheet,  the  mapping  is  applied  to  only  those  poles, 
and  the  poles  are  represented  alone  on  a  plane.  This  plane  really 
only  corresponds  to  the  one  sheet,  but  such  a  representation  removes 
confusion  concerning  other  poles  of  no  current  interest.  Hence,  this 
is  done  in  this  thesis.  An  example  is  the  mapping  in  Figure  1. 


The  systems  in  the  example  problems  which  follow  are  of 
order  one  half  and,  thus,  require  surfaces  of  two  sheets.  One  sheet 
is  taken  as 

50  =  {  z  |  -  k  <  arg(z)  <  k  } 

In  order  to  avoid  the  confusion  of  jump  discontinuities  experienced 
by  poles  moving  across  the  branch  cut  joining  cuts  defined  at  the 
extremes  of  the  arg(z)  range,  e.g.  0  and  4k  ,  this  thesis  uses  two 
half-sheets: 

S-i  =  {  z  |  -  2tc  <  arg(z)  <  -  n  } 

51  =  {  z  |  k  <  arg(z)  <  2k  } 

S-i  and  Si  are  joined  at  arg(z)=2n  and  arg(z)=  -2k.  S-i  is  joined  to 
So  at  the  branch  cut  where  arg(z)  =  -k,  and  Si  is  joined  to  So  at  the 
branch  cut  where  arg(z)  =  k  .  The  advantage  of  this  choice  of  sheets 
is  that  mapping  to  the  s1 -plane  is  done  by  halving  the  argument  of 
the  quantity  being  mapped  and  no  poles  move  across  the  branch  cut 
joining  the  extreme  argument  values.  Consequently,  the  jump 
discontinuities  in  the  argument  are  avoided.  In  summary,  the  poles 
lying  in  the  s1/n-p!ane  outside  the  sectors  corresponding  to  the  So 
Riemann  sheet,  considered  the  s1 -plane  in  this  thesis,  reside  on 
other  Riemann  sheets.  These  poles  have  been  shown  to  contribute 
terms  of  negative  powers  of  time  to  the  response  (1:92-104),  and 
these  terms  appear  to  correspond  to  creep  or  relaxation  response  in 
viscoelastic  materials  (2:748,  20:681).  As  such,  this  response  is 
stable  and  non-oscillatory.  Therefore,  the  poles  lying  on  the  So 
sheet  are  solely  responsible  for  the  system  stability  and  oscillatory 
motion,  because  their  response  is  exponential,  Ae®*,  which  can  be 
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oscillatory  and  stable  or  unstable,  accordingly  as  Re  co  £  0  or 
Re  co  >  0,  respectively. 


Finally,  to  complete  the  discussion  of  the  solution  for  system 

response,  the  method  of  actually  calculating  the  total  response  is 

necessary.  A  detailed  description  of  the  calculation  of  the  inverse 

— * 

Laplace  transform  of  the  state  vector  transform,  X(s),  is  given  at 
the  end  of  Example  1  in  Chapter  IV  after  the  reader  has  become  more 
familiar  with  the  pole  structure  of  systems  of  order  1/n.  The 
calculation  involves  the  evaluation  of  a  contour  integral  in  the 
s1-plane. 


1.  Formulate  the  system  of  order  1/n  (as  illustrated  in  the 
examples). 


2.  Determine  the  optimal  gain  matrix  K  through  an  algebraic  Riccati 
equation  solver,  such  as  RICCATI  in  MATRIXx  ,  by  using  the 
equivalent  first-order  system,  Eqs  (37-40). 

3.  Solve  the  eigenvalue  problem  for  the  eigenvalues,  X,  and 
eigenvectors,  ,  .in  the  s1/n-plane. 

4.  Map  the  poles,  X,  to  the  s1 -plane  poles,  co. 

5.  Calculate  the  total  response  (as  illustrated  in  Example  1). 


36 


1Y.  Example  Problems 


In  this  chapter  the  theory  is  demonstrated  by  numerical 
solution  by  finite  element  analysis  of  two  examples  of  simple 
controlled  structures.  An  iterative  solution  of  the  open  loop 
equation  of  motion  for  the  first  mode  eigenvectors  and  eigenvalues 
allows  cross-checking  the  open  loop  eigenvalues  and  eigenvectors 
computed  by  matrix  manipulation  with  MATRIXx-  The  first  example 
is  a  beam  of  viscoelastic  material,  pinned  to  allow  only  rotation  at 
both  ends.  Both  rotations  are  controlled.  The  second  example  is  a 
long,  thin  aluminum  rod,  clamped  at  one  end  and  damped  by  a 
viscoelastic  pad  at  the  other.  Displacement  is  constrained  to 
longitudinal  motion  and  is  controlled  at  the  damped  end. 

These  problems  illustrate  the  method  of  applying  the  control 
theory  to  structures,  as  well  as  the  technique  of  expanding  the 
structure  equation  of  motion  into  a  differential  system  of  order  1/n. 
The  first  problem  has  a  lot  of  symmetry  (damping,  control,  boundary 
conditions,  geometry),  while  the  second  does  not.  The  second  shows 
how  to  incorporate  both  viscoelastic  and  elastic  materials  into  a 
problem.  The  basic  concept  of  finite  element  analysis  does  not 
change;  viscoelastic  materials  are  included  in  the  element  equation 
of  motion  and  boundary  conditions.  The  problems  also  show  the 
technique  of  using  MATRIXx  to  solve  the  Riccati  and  eigenvalue 
equations. 

Example  Problem  1:  Simply-Supported  Beam 
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Problem  Statement:  Given  a  simply-supported  beam  (pinned  c*t 
both  ends  to  allow  only  rotation)  made  of  butyl  B252  rubber  as 
shown  beiow,  determine  the  open  and  closed  poles  to  illustrate  the 
effect  of  the  linear  feedback  theory  on  system  response,  as 
described  by  rotational  displacements  0i  and  02  at  nodes  1  and  2. 


b 

-H 
C 

X-section 


Figure  2.  Example  1.  Simply-supported  Beam 


Geometry: 

L  =  1.0  m 
b  =  0.1  m 
h  =  0.1  m 

I  =  bh3/12  =  8.33E-6  m4  (Moment  of  Inertia) 

E  =  a  time-dependent  modulus  to  be  described  later  in  the  Laplace 
domain 

Boundary  Conditions: 

vi  =  V2  =  0  (transverse  linear  displacements  of  nodes) 
ui  -  control  torque,  determined  by  optimal  control,  u*,  on  0i 


-H 


U2  =  control  torque,  determined  by  optimal  control,  u\  on  02 
Stress-Strain  Model 

As  Bagley  (3:920)  demonstrated,  Butyl  B252  rubber's  tensile 
stress-strain  relation,  can  be  written  in  terms  of  a  fractional 
derivative  model  of  order  a,  which  is  determined  by  experiment 
involving  sinusoidal  excitation  over  a  range  of  frequency  0  <co<104 
cps: 

otD+bdJort-EoeW+^Cfcd)  (2) 

where  o(t)  is  the  normal  stress,  e(t)  is  the  normal  strain,  a  =  m/n  (m, 
n  are  integers),  and  Eo  and  Ei  are  positive,  constants.  For  this  work 
b  equals  zero.  In  the  time  domain  an  operator,  PE,  acting  on  the 
strain,  e(t),  can  be  defined  such  that  Eq  (2)  can  be  written  as 

a(t)  =  PE(  e(t)  )  (45) 

In  the  Laplace  domain,  this  operator  can  be  described  as  a  funtion, 
E(s),  where  the  overbar  denotes  the  Laplace  transform  domain,  of  the 
Laplace  variable,  s.  This  can  easily  be  seen  if  one  takes  the  Laplace 
transform  of  the  general  stress-strain  model,  which  can  be  done 
since  the  model  is  linear.  Assuming  zero  strain  at  t=0, 

a(s)  =  Eo  e(s)  +-^-[  si(s)] 
s1  *“ 

a(s)  =  (E0  +  E^)  e(s) 

a(s)  -  E(s)  e(s)  (46) 


The  function  E(s)  is  not  a  transform;  it  corresponds  to  an  operator  in 
the  time  domain  and  not  to  a  function.  This  form  of  the  time- 


dependent  modulus,  i.e.  E(s),  is  useful  for  applying  the 
viscoelasticity  correspondence  principle  (5:42).  This  principle 
states  tnat  in  tne  Laplace  domain  the  time-dependent  modulus  can 
be  substituted  for  a  presumed-elastic  (constant)  modulus,  E,  in  an 
equation.  Because  the  operator  PE(e(t))  can  be  expressed  merely  as 
an  polynomial  expression  in  s  in  the  Laplace  domain,  it  can  be 
grouped  as  one  term,  as  shown  above.  The  correspondence  principle 
says  that  E  can  be  replaced  by  E(s)  in  the  Laplace  domain  when 
viscoelastic  loads  are  calculated  using  convolution  which 
corresponds  to  multiplication  in  the  Laplace  domain.  When 
viscoelastic  materials  are  used  in  finite  element  models,  the 
equation  of  motion  (EOM)  for  an  element  is  first  derived  assuming  a 
constant  modulus,  E.  The  correspondence  principle  is  then  applied  to 
the  transform  of  the  EOM.  The  resulting  equation  is  solved  'or  the 
response. 

Based  on  Madigowski  data,  Bagley  (3:920)  has  determined  a  and 
positive  shear  coefficients  Go  and  Gi  in  the  shear  stress-strain 
relationship  for  butyl  rubber  and  several  other  materials: 

x(t)  =  Goy(t)  +GiD?y(t)  (47) 

where  x  is  shear  stress  and  Y  is  the  shear  strain.  To  determine  Eo 
and  Ei  one  applies  the  basic  relation  between  Ey,  Young's  modulus, 
and  G,  shear  modulus,  i.e.  Ey  =  2(1+v)G,  where  v  is  the  material's 
Poisson's  ratio,  and  uses  Go  and  Gi  to  determine  Eo  and  Ei , 
respectively.  For  this  work  the  rubber  was  assumed  incompressible, 
so  v  =  0.5.  Therefore, 
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Butvl  B252  Rubber  Values 
p  =  0.92  g/cm3  =  920  kg/m3 
a  =  0.5 

G0  =  760,000  N/m2 
Gi  «  295,000  Ns1/2/m2 
E0  =  2,280,000  N/m2 
Ei  =  885,000  Ns1/2/m2 

These  figures  are  valid  for  0  <  frequency  <  104  Hz. 

Elastic  Formulation  of  Equation  of  Motion 

In  finite  element  modelling  the  general  equation  of  motion  for 
a  sturcture  is  formed  by  assembling  the  equations  of  motion  of 
individual  elements  into  one  large  system  equation.  The 
correspondence  principle  is  applied  to  the  element  equation  of 
motion  for  those  elements  containing  the  viscoelastic  material. 
These  modified  equations  are  assembled  into  the  system  equation  of 
motion,  and  the  EOM  is  solved  for  the  system  response. 

The  general  form  of  the  equation  of  motion  for  an  elastic 
finite  element  is 

md  +  kd  =  7  (48) 

where 

m  =  element  mass  matr  ix  (lumped  or  consistent) 
d  =  element  displacement  vector 
k  =  element  stiffness  matr  ix 
7  =  element  applied  nodal  loads 


Cook  (7:87)  gives  the  stiffness  matrix  for  an  elastic  beam  element 


as 

12  6L  -12  6L 
~  =  Ej_  6L  4L2  -6L  2L2 
L3  -12  -6L  1 2  6L 
.  6L  2L2  -6L  4L2  . 

where  E  is  the  presumed-elastic  modulus  and  the  displacement 
vector  is  arranged  by  alternating  translational  and  rotational 
displacements  for  each  degree  of  freedom,  i.e. 


The  applied  element  force  vector,  r,  brings  in  the  control  torques: 


Since  rotational  displacements  are  of  interest,  Cook's  consistent 


mass  matrix  must  be 

used: 

156 

22L 

54 

-13L 

M  =  m 

22L 

4L2 

13L 

-3L2 

420 

54 

13L 

156 

-22L 

.  -13L 

CM 

—1 

CO 

1 

-22L 

4L2 

where 


m  =  pAL  (A  =  cross-sectional  area) 
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Since  there  is  only  one  element  in  the  structure,  the  assembled 
matrices,  hence  the  equation  of  motion,  are  the  element  mass, 
stiffness,  and  force  matrices.  Applying  the  aisplacement  boundary 
conditions  reduces  the  system  to  two  degrees  of  freedom: 


pAL 

420 


or 


4L/  -3U 

-3L2  4L2 


ImoI  ei 
1  02(1)  I  L3 


M  j  ®l(t)  )  +  E  K*  j  0l(t) 
1  02(t)  )  I  02(t) 


4L2  2L2 
2L2  4L2 


Ui(t)  \ 

u2(t)  I 


I  01  (t) 

\  0  2  ( t) 


Ul(t)  | 

U2(t)  i 


(49) 


where 


m=?al' 


420 


4 
L  -3 


-3 

4 


£  =  2- 


2  1 
LI  2 


Because  E  will  be  replaced  by  E(s),  the  stiffness  matrix,  K,  has  been 
identified  as  EK'. 

Application  qf  lh£  Correspondence  Principle 

The  Laplace  transform,  denoted  by  capital  letters  of  functions, 
of  the  system  is 


Ms20(s)  +  EK'©(s)  =  U(s) 
where  s  is  the  Laplace  variable  and 


(50) 
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> 


0(S)  = 

I 


I  ©l(s) 
1  ©2<s) 


U(s)  = 


|U,(s)\ 

lUa(s)  j 


E(s)  replaces  E  and  gives  the  equation  of  motion  in  the  Laplace 
domain: 


[  Ms2  +  EiS1  72  K’  +  EqK'  0(s)  =  U(s)  (EOM) 


(51) 


Open  Loop  Response 

Recall  from  Chapter  III  that  the  Laplace  variable  s  corresponds 
to  the  eigenvalue  for  a  first-order  system.  Because  the  fractional 
derivative  transforms  have  s  to  fractional  powers,  an  eigenvalue  X 
is  defined  in  the  s1/2-plane  by  X  =  s1/2  ,  in  order  to  have  a 
polynomial  with  only  integer  powers.  The  eigenvector,  $,  is  given  by 
the  specific  value  of  ©( s)  for  the  eigenvalue  X.  Then  the  open  loop 
or  homogeneous  equation  of  motion  is 


MX  +  E-|AK'  +  E0K'  j  <t>  =  0  (52) 

which  can  be  solved  recursively  for  its  first  mode  eigenvalue,  Xi , 

— • 

and  eigenvector,  <>,  in  the  s1/2-plane  via  an  iterative  routine  (3:922- 
923): 


Xi(n  }Ei  K’  +  E0K' 


-1~  ~(n) 

M  di  = 


1 


4.( n  +1) 

'1 


-(n  +1) 
01 


(53) 


where  (n  )  and  (n  + 1)  are  iteration  steps.  The  routine  steps  are 
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1.  Estimate  A./  ,  <j>i  Quality  of  guess  only  affects  speed 

of  convergence.  For  this  thesis  1+i  was  simply  used  for  A.i(1). 

-d  ) 

Undamped  open  loop  first  mode  eigenvectors  were  used  for  4>i 

2.  Iterate  on  4>i  holding  A.i(n)  constant  until  the  eigenvector 

converges,  i.e.  changes  less  than  a  tolerable  difference. 

in  +-n  -’(n  +1)  ln  x  -*(n  ) 

3.  Use  new  Ai  ;  and  4>i  as  estimates  for  A./  ',  <j>i 

4.  Repeat  steps  2  to  3  with  until  A.i(n  }  converges. 

Because  it  is  a  quartic  equation,  there  are  four  values  of  ?Li(n+1) 
generated  each  iteration:  A1/n  =  A1/n{  cos  [(Q  +  2nk)/n]  +  i  sin  [(Q  + 
2itk)/n]  }  where  Q.  is  the  argument  of  A,  A  is  the  magnitude  of  A,  and  k 
is  an  integer  ranging  from  0  to  n-1.  One  value  must  be  chosen  for 
the  iteration.  Since  each  value  corresponds  to  the  possible  solution 
on  its  branch  or  sector  of  the  plane,  the  same  solution  must  be 
pursued  to  gain  convergence  to  the  valid  eigenvalue.  Identification 
of  branches  is  not  difficult  for  open  loop  poles  if  one  assumes  that 
the  pole  argument  does  not  vary  greatly.  Because  fourth-roots  are 
separated  by  ic/2  radians,  two  can  never  lie  in  the  same  quadrant. 
Therefore,  each  branch  can  be  identified  by  the  signs  of  the  real  and 
imaginary  parts.  To  maintain  the  same  branch,  always  iterate  with 
the  root  with  the  specified  signs.  This  works  well  for  open  loop 
poles,  but  closed  loop  poles  may  move  drastically  from  the  open  loop 
positions  and  into  other  quadrants.  Improved  algorithms  for  this 
situation  have  been  studied  currently  by  Captain  Michele  Devereaux 
(8). 
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Furthermore,  because  it  is  a  quartic  equation,  there  are  four 
distinct  eigenvalues  for  each  eigenvector.  The  algorithm  converges 
to  only  one  at  a  time,  so  the  iteration  must  be  completed  for  each. 
Eigenvalues  will  always  be  real  or  occur  in  complex  conjugate  pairs 
because  the  coefficients  are  real.  Moreover,  there  should  be  an  open 
loop  pole  in  each  quadrant,  because  the  viscoelastic  model  yields 
both  exponential  and  power  time  response.  The  right  two  quadrants, 
because  they  map  to  the  s1 -plane,  give  exponential  response  (stable 
or  unstable);  while  the  left  two  quadrants  of  the  s1/2-plane 
correspond  to  other  planes  in  which  poles  give  power  response. 
Hence,  one  set  of  poles  should  be  in  the  first  and  fourth  quadrants; 
and  one  set  should  be  in  the  second  and  third  quadrant.  There  may  be 
two  negative  real  poles  which  correspond  to  lying  on  planes  not 
shown  in  s1/2-plane  presentation.  This  procedure  must  be  modified 
slightly  for  higher  eigenvalues  (3:922-923,1:76-78). 

Applying  the  problem  values  to  the  EOM  matrices  gives 


M  -  0.0219  kg-m2 


-3 
4  . 


EiK- 


14.75  N-m 


2  1 

.1  2. 


E0K  =  38.0  N-m 


2  1 

.1  2- 


The  mass  matrix  M  is  multiplied  by  rad/sec2  in  the  equation  of 
motion.  Thus,  the  final  units  are  N-m,  i.e.  torque  which  is  balanced 
in  the  equation  of  motion. 


The  controls  analysis  package  (code)  MATRIXx(14)  was  used  to 
perform  the  iterative  solution.  MATRIXx  has  a  feature  called  a 
command  file  (14:7-3  -  7-6)  which  allows  a  series  of  commands  to 
be  stored  in  an  executable  file.  It  is  basically  the  equivalent  of 
program.  A  MATRIXx  command  file  entitled,  4th-order  BEAM 
Eigenvalue  Solver  or  4BEAMEVS.EXC,  was  created  to  solve  the 
recursive  equation  for  the  beam  open  loop  eigenvalues  of  the  first 
mode.  The  user  supplies  the  assembled  mass  and  stiffness  matrices, 
as  well  as  estimates  for  the  first  mode  eigenvector  and  one  of  the 
eigenvalues.  The  program  is  interactive,  and  the  user  is  asked  to 
identify  the  desired  root  from  the  output  and  input  it  for  iteration. 
The  difference  between  past  and  present  step  values  of  the  desired 
eigenvalue  is  displayed  for  both  real  and  imaginary  parts,  and  the 
user  is  queried  about  continuing  the  iteration  on  the  branch.  This 
could  be  automated,  if  desired.  For  this  work  the  iteration  was 
continued  until  the  difference  between  past  and  present  step  values 
of  the  desired  eigenvalue  was  less  than  10'5on  both  real  and 
imaginary  parts. 

Because  of  the  symmetry  of  the  problem,  especially  the 
damping,  the  damping  is  expected  to  be  proportional  damping,  i.e.  the 
damping  matrix  EiK  can  be  written  as  a  linear  combination  of  the 
mass  and  stiffness  matrices  (19:197).  By  inspection  it  is  obvious 
that  EiK  can  be  written  in  terms  of  the  EoK.  In  this  case,  the 
equations  comprising  the  system  can  be  decoupled,  and  the 
eigenvectors  are  the  same  as  those  for  elastic  system.  The  elastic 
first  mode  (19:213)  is  known  to  be 
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and  the  first  guess  of  the  eigenvalue  was  Xi  *  1  +  i.  Two  conjugate 
pairs  for  the  first  mode  result  from  the  iteration: 

Jti  =  2.9747  ±  4.1 1 50i  and  -2.9747  ±  0.8740i 


The  poles  with  the  positive  real  part  correspond  to  vibrational 
motion,  while  the  poles  with  the  negative  real  part  correspond  to 
relaxational  motion. 

The  iteration  also  gave  the  elastic  first  mode  eigenvector, 

Eq  (54),  for  all  eigenvalues.  This  occurs  because  of  the  proportional 
damping. 

These  values  can  be  checked  because  the  first  eigenvector  is 
known.  By  premultiplying  by  the  transpose  of  the  jth  eigenvector 
Eq  (52)  reduces  to  a  fourth-order  scalar  polynomial  in  Xj: 


-T 


X?M 


+  Xj  E-|K+  E0K 


♦j 


=  0 


By  using  a  root  finding  routine,  such  as  the  ROOTS  command 
(14:4-16)  in  MATRIXx  ,  the  eigenvalues  are  found  to  be  the  same.  The 
elastic  second  mode  is  also  known: 


and  the  eigenvalues  can  be  found  by  the  same  method: 


X2  =  7.1139  1  1 1 ,0272i  and  -11.6288,-2.5989 
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The  poles  with  the  positive  real  part  correspond  to  vibrational 
motion,  while  the  poles  with  the  negative  real  part  correspond  to 
relaxational  motion. 


The  general  problem  can  be  posed  in  a  fractional  state-space 
by  a  1/2th-order  system  (3:921-922).  A  fractional  state-space  is 
defined  as  a  state-space  in  which  fractional  as  well  as  integer 
derivatives  of  variables  are  included  in  the  state  vector,  e.g. 

3  r>~* 

D,  0(t) 

o  ( t\  _  .  0  ( t) 

d!5(«' 

0(t) 
or 

s3  ^©(s) 
s1 0(s) 
s1  *0(s) 

0(s) 

whe4m  states: 

‘  0  0  0  M  1  s3C®<s>  0  0  -M  0  l(s320(s) 

sw>  0  0  M  0  s1  0(s)  +  0  -M  0  0  s1  0(s) 

0  M  0  0  s1  e  0(s)  -M  0  0  0  s120(s) 

-  M  0  0  EiK  J  ©(S)  L  0  0  0  E0K  J  0(S) 
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U(s) 


0 

0 

0 


where  0  = 

'o  0 ' 

and  1  = 

1  0  ' 

.0  0  . 

.0  1  . 

This  is  rearranged  to 


s1  *0(s)  =  A  0(s)  +BU(s) 


where 


0 

0 

0  M 

0 

0 

-M 

0 

A  = 

• 

0 

0 

M  0 

0 

-M 

0 

0 

0 

M 

0  0 

-M 

0 

0 

0 

M 

0 

0  EiK  _ 

- 

0 

0 

0 

E0K  . 

- 

0 

0 

0  M 

0 

B  = 

0 

0 

M  0 

0 

0 

M 

0  0 

0 

M 

0 

0  EtK  . 

.  1 

_ 

(56) 


(57) 


To  complete  the  system,  it  is  assumed  that  all  states  are  outputs 
and  that  there  is  no  feedforward  to  the  output,  i.e. 

Y(s)  =  C0(s)  +  D  U(s) 
where 

C=  4m  x  4m  identity  matr  ix 

D  ■  m  x  q  null  (zero)  matr  ix  (q  =  number  of  contro  I  inputs) 

For  the  specific  problem  m=2  and  q  =  2,  so  the  system  has 
eight  states  and  two  control  inputs. 
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Open  Loop  Poles  via  MATRIXy 


If  the  centrel  vector  is  set  tc  zerc,  i.e.  0,  in  the  1/nth-crder 
system,  Eq  (57)  reduces  to  an  eigenvalue  problem.  However,  there 
will  be  m  sets  of  eigenvalues,  where  m  is  the  number  of  degrees  of 
freedom  again.  Eigenvalues  and  eigenvectors  of  the  A  matrix  were 
determined  by  the  EIG  command  of  MATRIXx  (14:4-23).  These 
corresponded  to  the  first  two  modes.  The  first  mode  eigenvalues 
and  eigenvectors  determined  by  the  recursive  solution  of  the  open 
loop  equation  of  motion,  Eq  (52),  match  those  determined  by 
MATRIXx. 

Closed  Loop  Poles  via  MATRIXy 

Solution  for  closed  loop  poles  and  eigenvectors  is  quite  easy 
on  MATRIXx.  First,  the  system  must  be  posed  for  the  Riccati 
Equation  solver  in  the  appropriate  An  -An  B  form,  as  indicated  in 
Chapter  III: 

a(t)  =  An&(t)  +  An'1Bu(t ) 

Therefore,  for  this  example 

S(t)  =  A2£(t)  +  ABu(t)  (58) 

The  MATRIXx  Riccati  solving  command,  called  RICCATI 
(14:14-19),  was  used  to  determine  the  optimal  gain  matrix,  P,  such 
that 
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u*  ( t)  =  -R^V'1  Pa(t)  (59) 

However,  the  routine  also  output  the  overall  gain  matrix,  KC,  where 

KC  =  r'bV^P 

and  this  matrix  was  used  in  the  calculation  of  closed  loop  poles. 
RICCATI(S,R3NS)  requires  three  inputs:  a  system  matrix  (S),  a 
performance  index  weighting  matrix  (RQ,  and  the  number  of  states 
(NS)  in  the  state  vector.  These  are 

~n  Tn* 1  ~ 

S  _  A  A  B 

.  C  D  . 

RQ  =  5 

.  Nq  R 

NS  =  4m 

where 

Nq  =  4m  x  4m  null  matr  ix  (same  dimensions  as  Q) 

Nr  =  4m  x  q  null  matr  ix  (same  dimensions  as  R) 
m=  2  for  this  problem 
q  =  2  for  this  problem 

For  this  problem  m=  2  and  q  *  2.  The  closed  loop  system,  Eq  (58), 
becomes  an  eigenvalue  problem  once  the  gain  matrix  for  the  optimal 
steady-state  linear  feedback  is  known.  The  Laplace  transform  of 
Eq  (58)  is 

s1  *§(s)  =  (  A  -BR'1  BTA  Tp)q(s)  (59) 
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For  this  specific  problem  the  inputs  to  RICCATI  are 


~  A2  AB 
S  = 

.  C  D  . 

0  0  0  0  0 
_  0  10  0  0 
RQ  =  ooooo 
0  0  0  10 
io  oo  r. 

NS  ■  8 


MATRIXx  then  gives  P  and  KC: 


p  = 

1000  * 

0.0000 

o.ooco 

•0.0001 

0.0001 

0.0002 

-0.0001 

0.0025 

-0.0008 

0.0000 

0.0000 

-0.0001 

0.0001 

0.0002 

-0.0001 

0.0025 

-0.0008 

0.0000 

0.0000 

0.0001 

-0.0001 

-0.0001 

0.0002 

-0.0008 

0.0025 

-0.0001 

0.0001 

0.0005 

-0.0003 

-0.0025 

0.0008 

0.0000 

0.0000 

0.0001 

-0.0001 

-0.0003 

0.0005 

0.0008 

-0.0025 

0.0000 

0.0000 

0.0002 

-0.0001 

-0.0025 

0.0008 

0.0186 

0.0021 

0.0409 

-0.0091 

-0.0001 

0.0002 

0.0008 

-0.0025 

0.0021 

0.0186 

-0.0091 

-0.0409 

0.0025 

-0.0008 

0.0000 

0.0000 

-0.0409 

-0.0091 

2.2763 

1.7681 

-0.0008 

0.0025 

0.0000 

0.0000 

-0.0091 

-0.0409 

1.7681 

2.2763J 

KC  = 

_ 

-0.6380 

0.1986 

5.6916 

0.6239 

48.2600 

-26.8662 

0.0088 

-0.0044 

0.1986 

-0.6380 

0.6239 

5.6916 

-26.8662 

48.2600 

-0.0044 

0.0088J 

MATRIXx's  EIG  command  (14:4-23)  solves  the  closed  loop  eigenvalue 
problem  for  the  closed  loop  poles  and  eigenvectors  of  both  modes. 
The  poles  are  given  in  a  jumbled  order: 

7.2744  +11.2328i 
7.2744  -11.2328i 


-0.2318  +  3.0553i 
2.9596  +  4.2008i 
2.7549  +  4.6337i 
-0.2318  -  3.0553i 
2.9596  -  4.2008i 
2.7549  -  4.6337i 

The  command  EIG  does  not  identify  which  closed  poles  correspond  to 
which  open  loop  poles.  In  general,  if  the  open  loop  poles  are  known 
for  each  mode,  the  closed  loop  poles  for  each  mode  can  be  identified. 
Starting  from  the  open  loop  system,  feedback  gain  is  increased  from 
zero  to  full  value  in  small  enough  steps  to  see  the  poles  move  from 
open  loop  to  closed  loop  positions.  The  tracks  map  each  open  loop 
pole  to  its  closed  loop  position.  In  controls  analysis  this  plot  is 
ca! led  a  rcct-lccus.  The  eigenvalues  for  each  mode  of  the  specific 
problem  are 


Xi  =  2.9596  ±  4.2008i  and  -0.2318  ±  3.0553i 
\2  =  7.2744  ±  1 1 .2328i  and  2.7549  ±  4.6337i 


The  eigenvectors,  as  determined  by  the  EIG  command,  remain 
the  same  as  the  open  loop  eigen  vectors,  which  are  the  elastic 
eigenvectors  Eqs  (54)  and  (55).  This  makes  identification  of  the 
poles  easy. 

The  closed  loop  poles  can  be  determined  by  creating  a  scalar 
equation  using  the  known  eigenvectors,  as  in  the  open  loop  case.  The 
control  vector  is  replaced  in  the  system,  Eq  (58),  by  the  optimal 
control  vector,  which  is  KC  <}> .  The  vector  is  brought  to  the  left- 
hand  side  of  the  equation.  The  scalar  equation  of  motion  is  given  by 
the  bottom  m  rows,  where  m  is  the  number  of  degrees  of  freedom, 
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of  the  system  equation  and  now  has  non-zero  coefficients  for  X3  and 
X2  and  becomes 


— T 

*j 


XjM 


+  XiKCi, 


4-  X.?KCor 


.  1  I 


O  .  IS 


\  I 


where  KCqm  is  the  qth  set  of  m  columns  of  KC.  The  first  and  second 
mode  eigenvalues  are  found  to  be  the  same  as  the  quoted  above.  The 
Bagley  and  Torvik  algorithm  can  be  modified  to  solve  this  new 
eigenvalue  problem.  Unless  the  eigenvectors  are  known,  the  only 
other  means  of  determining  which  eigenvalues  belong  to  which 
modes  is  the  iterative  method  of  increasing  gain. 


Measures  of  Change  in  Stability 

A  quick  check  for  increased  stability  can  be  done  by  comparing 
pole  positions:  are  closed  loop  poles  further  left  (i.e.  real  part)  than 
open  loop  poles,  and  do  closed  loop  poles  have  a  larger  ratio  of 
imaginary  to  real  components? 

In  the  s1  -plane  the  real  component  equals  -Cwn,  where  t,  is  the 
damping  ratio  and  o)n  is  the  natural  frequency.  The  amplitude  of  the 
response  is  driven  by  e-C«nt.  Consequently,  pole  movement  further 
left  on  the  real  axis  results  in  a  more  negative  exponent:  and,  the 
system  amplitude  dies  faster.  A  measure  of  this  motion  is  given  by 
the  difference  between  open  and  closed  loop  real  components. 
Mapping  poles  between  the  s1/2-  and  si-planes  will  effect  the 
magnitude  of  the  value  but  will  not  change  the  sign  which  identifies 
the  d'rection  of  motion  (plus  means  left  and  negative  right). 
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The  imaginary  component  of  the  pole  in  the  s1 -plane  is  cod,  the 
damped  frequency.  Therefore,  the  ratio  of  imaginary  to  real 
components  gives  a  measure  of  both  the  damping  ratio  and  the 
damped  frequency  with  respect  to  the  natural  frequency.  Moreover, 
the  cosine  of  angle  between  the  negative  real  axis  and  a  line  drawn 
between  the  pole  and  the  origin  is  the  damping  ratio,  The  larger 


the  angle  the  smaller  is  the  damping  ratio,  which  means  less 
damping.  The  ratio  of  real  to  imaginary  parts,  thus,  is  also  a 
measure  of  £.  Since  stable  response  poles  lie  in  the  two  s1/2-plane 
sectors  rc/4  <  9  <  rc/2  and  -rc/2  <  0  <  -rc/4  ,  a  higher  ratio  puts  a 
pole  closer  to  the  imaginary  axis  (larger  tan  0;  hence,  larger  0)  than 
another.  Therefore,  it  maps  closer  to  the  real  axis  in  the  s1 -plane 
and  has  a  lower  damped  frequency  and  higher  damping  ratio, 

Because  the  poles  in  the  left  half  s1/2-plane  do  not  lie  on  the 
s1  -plane  and,  thus,  do  not  contribute  to  typical  stability  measures, 
their  movement  with  respect  to  stability  is  less  clear.  These  poles 
cause  negative-power  time  response  which  decays  slower  than  an 
exponential  response.  If  they  move  right  or  increase  their 
imaginary-to-real  ratio,  they  move  toward  the  s1  -plane.  If  they 
move  onto  the  left  half  of  the  s1 -plane,  their  response  is  more 
stable  than  before,  because  the  poles  now  cause  decaying 
exponential  response.  Therefore,  motion  right  and  increasing  ratios 
of  poles  in  the  left  half  s1/2-plane  would  seem  to  indicate  increased 
system  stability.  Under  these  rules  these  two  measures,  difference 
between  real  components  and  imaginary-to-real  ratio,  give 
measures  of  the  change  in  stability  of  the  system. 
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Comparison  of  Open  and  Closed  Loop  Poles 

The  open  loop  eigenvalues  of  the  first  and  second  modes  are  shown 
in  Figure  3. 
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Figure  3.  Example  1  Open  Loop  Poles  for  Both  Vibrational  and 
Relaxation  Modes  1  and  2. 


and  the  closed  loop  eigenvalues 


are  shown  graphically  in  Figure  4. 
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Figure  4.  Example  1  Closed  Loop  Poles  for  Both  Vibrational  and 
Relaxation  Modes  1  and  2. 


Real  components  of  the  open  and  closed  loop  eigenvalues,  Re(O) 
and  Re(C)  respectively,  the  difference  between  the  open  and  closed 
loop  real  parts,  denoted  by  Re(0^-Re(C),  and  ratios  of  the  imaginary 
to  real  parts,  denoted  by  l/R,  for  the  first  two  modes  are  tabulated 
in  Table  I.  A  negative  number  in  the  Re(0)-Re(C)  column  indicates 
the  poles  have  moved  right. 


Table  I.  Stability  Measures  for  Example  1 


Re(C) 

Re(0)-Re(C) 

Open  l/R 

Closed  l/R 

1 

2.9747 

2.9596 

0.0151 

1 .3833 

1 .4194 

-2.9747 

-0.2318 

-2.7429 

0.2938 

13.1808 

2 

7.1  139 

7.2744 

-0.1605 

1 .5501 

1  .5442 

-1  1  .6288 

2.7549 

-14.3837 

0.0 

1  .6820 

-2.5989 

2.7549 

-5.3538 

0.0 

1.6820 
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The  data  from  this  table  shows  that  optimal  control  does  not 
necessarily  mean  most  stable  for  all  poles.  The  first  two 
eigenvalues,  which  represent  the  exponential  response  of  the  s1- 
plane,  of  the  first  mode  increased  in  stability,  as  indicated  by  the 
motion  left  of  the  real  part  (positive  difference  Re(0)-Re(C))  and  by 
the  increase  of  the  imaginary-to-real  ratio,  l/R.  However,  the 
positive-real-part  eigenvalues  of  the  second  mode  destabilized 
slightly,  as  shown  by  the  motion  right  of  the  real  part  (negative 
difference  Re(0)-Re(C))  and  a  decrease  in  the  l/R  ratio.  The  second 
mode  poles  which  have  negative  real  parts  moved  right  far  enough  to 
enter  the  region  corresponding  to  the  left  half  s1  -plane.  Therefore, 
these  poles,  which 

formerly  gave  non-oscillatory  relaxation  response,  now  give  a 
vibrational  response!  Furthermore,  the  poles  lie  nearly  on  top  of  the 
positive-real-part  poles  of  the  first  mode.  Because  of  the  first  and 
second  mode  eigenvectors,  Eqs  (54)  and  (55),  there  is  constructive 
interference  at  node  1  and  destructive  interference  at  node  2. 

Finally,  the  negative-real-part  eigenvalues  of  the  first  mode  also 
destabilized  by  moving  right  but  has  increased  its  damping  ratio  and 
lowered  its  damped  frequency,  as  indicated  by  the  increased  l/R 
ratio. 


Most  of  the  closed  loop  poles  are  not  very  far  from  the  open 
loop  poles.  It  would  appear  that  these  values  are  already  fairly 
optimal  under  the  current  evenly-weighted  performance  index.  The 
problem  appears  to  be  heavily,  perhaps  overly,  damped,  and  the 
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destabilization  of  the  one  pole  may  represent  a  desire  for  the 
system  to  be  less  damped.  By  destabilizing  the  pole  the  system 
would  appear  more  flexible. 

The  system  is  stable  with  the  optimal  feedback,  which 
produces  significant  changes  in  the  response.  However,  this  system 
is  not  very  indicative  of  large,  flexible  space  structures  because  of 
the  geometry,  i.e.  this  system  is  heavy  weight.  Furthermore,  there 
is  so  much  damping  material  that  the  system  is  probably  overdamped 
and  does  not  need  much  feedback,  if  any.  Therefore,  a  problem  more 
representative  of  the  structures  of  interest  is  needed. 

Modified  Beam  Geometry  to  Represent  a  Flexible  Space  Structure 

Such  a  problem  can  be  achieved  by  a  modification  of  the 
geometry  to  represent  a  more  slender  beam.  Therefore,  let  b  =  h  = 
0.01m  .  This  divides  M  by  100  and  divides  K'  by  10*.  The  open  loop 
eigenvalues  are 

A.i  ,ol  =  0.8936  ±  1 .0333i 
-0.8936  ±  0.7276i 

and 

\2,oi  =  1.9573  ±  2.5320i 
-1.9573  ±  1.11 83i 


and  are  shown  graphically  in  Figure  5. 


♦  VIBRATIONAL  1 
o  VIBRATIONAL  2 

H 

3  ■  RELAXATION  1 

□  RELAXATION  2 


Figure  5.  Modified  Example  1  Open  Loop  Poles  for  Both  Vibrational 
and  Relaxation  Modes  1  and  2. 


By  following  the  same  procedures  outlined  earlier,  the  gain  matrix 
•  KC  is  determined  to  be 


KC  = 

-0.0007  0.0001  1.0009  -0.0007  -0.0036  -0.0014  0.9924  -0.0038 

0.0001  -0.0007  -0.0007  1.0009  -0.0014  -0.0036  -0.0038  0.9924J 


The  closed  loop  eigenvalues  are 

\i.cl  =  0.2583  ±  25.5363i 
0.0000  ±  I.OOOOi 


and 


A.2.CL  =  1.3870  ±  67.5522i 
0.0000  ±  I.OOOOi 


and  are  shown  graphically  in  Figure  6. 
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Figure  6.  Modified  Example  1  Closed  Loop  Poles  for  Both  Vibrational 
and  Relaxation  Modes  1  and  2. 


The  stability  measures  for  the  poles  are  shown  in  Table  II. 


Table  II.  Stability  Measures  for  Modified  Example  1 


m 

mm 

mil 

Open  l/R 

Closed  l/R 

1 

0.8936 

0.2583 

0.6353 

1.1563 

98.8630 

-0.8936 

0.0 

-2.7429 

0.8142 

oo 

2 

1  .9573 

1 .3370 

0.5703 

1 .2936 

48.7038 

-1 .9573 

0.0 

-1  .9573 

0.5713 

03 

The  data  in  Table  II  shows  stability  increases  with  the  optimal 
feedback  control  for  the  positive-real-part  poles  of  both  first  and 
second  modes.  These  poles  have  moved  left  (positive  Re(0)-Re(C)) 
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and  have  dramatically  increased  their  imaginary-to-real  (l/R) 
ratios.  The  latter  indicates  a  significant  increase  in  damping  ratio, 
£  ,  as  will  as  a  significant  decrease  in  the  damped  frequency.  Once 
again,  the  negative-real-part  poles  moved  right  (negative  Re(O)- 
Re(C))  as  well  as  increased  their  l/R  ratios.  In  fact,  these  poles 
achieved  the  maximum  value  of  the  l/R  ratio.  They  have  moved  to 
the  imaginary  axis  and  will  map  to  the  negative  real  axis  in  the  s1- 
plane.  Therefore,  these  poles,  which  previously  gave  relaxation 
(power-law)  response,  now  give  purely  monotonically  decaying 
exponential  response. 

The  positive-real-part  poles  produce  the  characteristic 
damping  and  frequency  parameters  for  the  oscillatory  motion.  The 
eigenvalues  in  the  s1/2-plane  and  s^plane,  denoted  by  EV  s1/2  and 
EV  s1  respectively,  and  oscillatory  motion  parameters,  which  are 
damping  ratio  (£),  natural  frequency  (ton),  and  damped  frequency  (co<j) 
,  are  shown  in  Table  III. 


Table  III.  Example  1  Oscillatory  Motion  Parameters 


EV  s1'2 

EV  si 

c 

®  n 

tod 

Open  Loop 

0.8936  ±  1.0333i 

-0.2692  ±  1.8467i 

0.1442 

1 .8669 

1  .8467 

1.9573  ±  2.5320i 

-2.5800  ±  9.91 18i 

0.2519 

10.2422 

9.91  18 

Closed  Loop 

0.2583  ±  25. 5363i 

-652.04  ±  1 3.1 92i 

0.9998 

652.17 

13.192 

1.3870  ±  67.5522i 

4561.4  ±  1 87.39i 

0.9992 

4565.1 

187.39 

The  data  shows  a  tremendous  improvement  in  system  response.  The 
damping  ratio,  C.  becomes  nearly  equal  to  one,  i.e.  critical  damping, 
which  causes  the  system  response  to  approach  equilibrium  faster 
than  any  other  damped  response.  The  resonant  frequencies  have  been 
driven  up  to  much  higher  frequencies.  Therefore,  the  oscillatory 
motion  dies  much  faster,  because  e-C^n1  governs  the  decay  of  the 
envelope  of  the  motion  amplitude.  Also,  the  ability  to  excite  the 
structure,  by  simple  handling  or  operation,  has  been  reduced,  e.g.  it 
would  not  be  excited  by  merely  walking  on  it.  Furthermore,  the 
resonant  mode  separation  has  been  increased  from  8.375  rad/sec  to 
3913  rad/sec.  This  effectively  isolates  one  mode  from  another,  and 
controllers  can  neutralize  one  mode  with  less  excitation  on  another. 
These  results  show  conclusively  that  the  optimal  linear  feedback 
control  utilizing  fractional  state-space  can  successfully  and 
significantly  increase  system  stability. 

Contribution.  oi  Each  Derivative  to  Pole  Structure 

Structures  have  been  controlled  by  position,  velocity,  and 
acceleration  feedback  in  the  past.  What  benefit  lies  in  also  feeding 
back  fractional  states?  Insight  to  the  answer  can  be  gained  by 
feeding  back  only  one  derivative  at  a  time. 

For  example,  feeding  back  only  velocity  (using  the  optimal 
gains  determined  for  the  fractional  state  model)  yields  poles  further 
left  than  those  from  feeding  back  the  entire  state  matrix: 

0.0007  ±  25.5571  i  and  0.0022  ±  67. 5737i  which  map  to 
-653.16  ±  0.0358i  and  -4566.2  ±  0.2973i.  Although  these  have  a 
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better  damping  ratio,  they  do  not  minimize  the  performance  index. 
More  highly  damped  is  not  necessarily  optimal.  However,  this  does 
show  that  the  velocity  is  a  major  player  in  shifting  the  poles  left. 

If  the  3/2-derivative  is  fed  back  alone,  the  s1/2-plane  poles 
are  1.0466  ±  1.0243i  and  2.9034  ±  2.3290i .  The  poles  have  moved 
right  and  have  even  become  unstable!  Therefore,  tne  3/2-derivative 
would  counteract  some  of  the  strong  stabilizing  influence  of  the 
first  derivative,  so  that  the  performance  index  is  minimized  further 
than  by  using  velocity  alone.  This  last  statement  derives  from  the 
fact  that  the  optimization  of  the  gain  did  not  zero  the  gain  for  the 
3/2-derivative.  If  the  poles  of  the  velocity-only  feedback 
minimized  the  performance  index,  then  this  zeroing  would  occur. 

If  both  3/2-derivative  and  first  derivative  are  fed  back,  then 
the  poles  become  nearly  the  final  optimal  values:  0.2590  ±  25.555 8i 
and  1.3892  ±  67.5595i  where  0.2583  ±  25.5363i  and 
1.3870  ±  67.5522i  are  the  optimal  poles  in  the  s1/2-plane.  Hence, 
these  derivatives  are  responsible  for  the  majority  of  the  movement 
of  the  poles. 

If  the  other  derivatives  are  fed  back  individually,  a  pattern 
emerges.  The  odd  half-derivatives  (1/2,  3/2,  etc)  move  the  poles 
right,  and  the  even  half-derivatives  (0,  1,  etc)  move  the  poles  left. 
Furthermore,  the  higher  the  derivative  the  greater  is  its  effect  in 
each  category.  So,  a  3/2-derivative  moves  the  poles  more  than  a 
1 /2-derivative.  The  gains  of  each  are  traded  off  to  minimize  the 
performance  index  with  the  lower  derivatives  acting  as  fine  tuners. 
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Because  the  higher  derivatives  result  in  higher  order 
polynomials  in  the  feedback  transfer  function,  they  can  contribute 
more  poles  to  move  the  root-locus  left  in  the  s1 -plane.  Take,  for 
example,  a  general  negative  feedback  system  with  plant  transfer 
function  G{s)  and  feedback  transfer  function  H(s)  in  the  s1/2-plane. 
The  characteristic  equation  is  1  +  G(s)H(s)  =  0.  If  h(t)  is  a 
1 /2-derivative,  then  its  transform  in  the  s1/2-plane  is  H(s)  =  As, 
where  A  is  a  constant.  H(s)  can  generate  only  one  pole  (root  of  the 
characteristic  equation)  for  the  system,  because  H(s)  increases  the 
orderof  the  characteristic  equation  by  only  one.  However,  if  h(t)  is  a 
3/2-derivative,  then  H(s)  equals  As3  and  generates  three  poles 
because  it  raises  the  order  by  three.  Location  of  the  poles  moves  the 
root-locus  left  or  right,  but  more  poles  will  cause  more  effect  than 
less.  Therefore,  the  higher  derivatives  have  more  of  an  effect  on  the 
pole  structure. 


Calculation  of  Total  Response 

Finally,  to  complete  the  discussion  of  system  response  and 
stability,  the  method  of  actually  calculating  the  total  response  is 
necessary.  This  will  hopefully  clarify  how  the  poles  which  have 
been  neglected  so  far  (i.e.  those  not  on  the  s1 -plane)  are  accounted 
for  in  the  total  response. 

Bagley,  in  his  dissertation  (1),  and  Bagley  and  Torvik  (2:747- 
748)  describe  the  calculation  of  a  system  response  to  impulsive 
loading.  Fundamentally,  the  inverse  Laplace  transform  of  the 
displacement  vector,  L-1[X(s)],  where  s  is  the  Laplace  variable  and 
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the  capital  X  indicates  the  Laplace  transform  of  the  state  vector, 
must  be  calculated.  By  definition,  the  inverse  transform  is  a 
contour  integration  in  the  s1 -plane,  as  shown  in  Figure  3.  Segments 
3,  4,  and  5  are  required  because  of  the  branch  cut  on  the  negative 
real  axis. 


+Re 


Figure  7.  Integration  Contour  to  Calculate  Response 


The  inverse  Laplace  transform  of  X(s)  is  defined  as 

n+i~ 

L  1  (X(s)]  =  — I  esl  X(s)  ds. 

2*'  Li. 

The  segment  of  the  contour  which  produces  the  inverse  transform 
L-1[X(s)j  is  segment  1.  The  transform  X(s)  can  be  written  as  a 
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summation  of  the  4m  eigenvalues  and  eigenvectors,  where  m  is  the 
number  of  degrees  of  freedom. 


v  i,TF(s)  (l-*s“)  ; 

la  - 7— - 

n=1  mn  (S1  ^Xn) 


where 

bn  =  nth  eigenvector  of  the  non-expanded  equation  of  motion 

F(s)  =  transform  of  the  non-expanded  applied  force  vector 

bn  =  nth  eigenvector  of  the  expanded  equation  of  motion 

M.  =  the  expanded  mass  matr  ix  which  is  the  coefficient  matrix 
associated  with  the  s1/2  term  in  Eq  (56),  an  early  form 
of  the  expanded  system,  Eq  (57) 


—  T  ~  - 

mn=bn  M_  bn  =  a  modal  constant 


Since  both  open  and  closed  systems  have  the  form  of  an  unforced 
(homogeneous)  system,  the  response  is  determined  by  examining  the 
impulse  response.  This  also  assumes  that  all  initial  conditions  are 
zero.  Therefore,  F(s)  =  1,  vector  of  ones,  because  the  impulse  is 
applied  to  all  nodes.  Also  recall  that  b=0  for  this  problem.  Another 
important  point  is  that  there  are  eight  eigenvalues  and  eigenvectors 
from  the  expanded  equation  of  motion  but  only  two  modes.  The 
eigenvectors  of  the  expanded  equation  of  motion  are  distinct,  but  the 
mode  shapes  of  physical  deflection  are  not.  If  the  eigenvectors  had 
complex  elements,  we  would  have  seen  complex  conjugate  pairs  of 
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eigenvectors.  This  is  shown  in  Bagley's  and  Torvik's  example  in 
their  paper  (2:746). 

By  including  the  transform  in  the  contour  integration  and  using 
the  residue  theorem  from  complex  calculus,  the  transform  can  be 
written  in  terms  of  a  sum  of  the  integrals  for  the  other  five 
segments  and  a  sum  of  the  residues: 

.  r1-  -  6  f 

-1—  est  X(s)  ds  =  — 1—  Y  est  X(s)  ds  +  Y  b, 

/r-i-  2 Kik~Jk  f  ' 

— * 

where  bj  are  the  residues  of  the  poles  enclosed  by  the  contour.  The 
residues  can  be  shown  to  produce  the  exponential  response  related 
to  the  s1 -plane  poles.  Bagley  (4:141-143)  has  shown  that  the 
integrals  for  segments  2,  4,  and  6  are  zero,  when  radii  of  2  and  6 
increase  toward  infinity  and  the  radius  of  4  approaches  zero.  He  has 
also  shown  that  over  long  time  the  integrals  for  segments  3  and  5 
asymptotically  approach  t-n  for  t  large,  where  ri  >  1.  It  is  in  these 
integrals  that  the  effect  of  the  poles  not  on  the  s1 -plane  becomes 
included  in  the  overall  response.  Only  here  are  they  included  in  the 
response  calculation  because  of  the  form  of  X(s).  As  has  been  noted 
earlier,  this  portion  of  the  response  correlates  to  the  power-law 
stress  relaxation  observed  by  Nutting  (11)  in  deformed  materials. 
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Figure  8.  Example  2.  Fixed-Shear  Rod  Diagram 


Geometry  and  Properties 
L=  1  m 

b  =  h  =  0.01 m 
p  =  2796  kg/m3 
E  =  6.894  x  ICO  N/m2 

Ac  =  0.0001  m  (area  of  contact  between  shear  pad  and  rod) 
th  =  0.01  m 

G(s)  =  Go  +  Gis1/2  (shear  modulus  defined  in  the  Laplace  domain) 

Go  =  7.6  x  105  N/m2 

Gi  =  2.95  x  105  Ns1/2/m2 

Boundary  Conditions  and  Applied  Forces 

di  -  0 

t  =  shear  stress  applied  by  pad  on  node  3 
u  =  control  force  applied  at  node  3 
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The  formulation  follows  like  any  other  finite  element 
formulation  again  with  the  added  time-dependent  modulus  as  in 
problem  1 .  The  basic  equation  of  motion  for  an  element  is  again 


md  +  kd  -  F 


(39) 


From  Cook  the  elastic  rod  element  mass  and  stiffness  matrices  are 


pAL [ 2  1 


k=££  1  -1 

L  L  -1  1 


The  derivation  of  the  general  equation  of  motion  by  summation 
of  forces  on  nodes  is  important  to  review  here,  in  order  to 
understand  the  inclusion  of  the  shear  force  boundary  condition  in  F. 
Consider  node  3,  shown  in  Figure  9. 


D'Alembert 


Node  3 


Shear 


Figure  9.  Free  Body  Diagram  of  Node  3  in  Example  2 
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A  positive  displacement  of  the  node  is  assumed  to  stretch  the  rod, 
like  a  spring,  so  the  force  on  the  node  is  in  the  negative  direction. 
Furthermore,  the  rod  mass  is  lumped  on  the  node,  so  there  is  inertia. 
This  creates  the  D’Alembert's  force  resisting  the  node  motion. 
Finally,  the  shear  force  created  by  the  shear  pad  is  assumed  to  be  a 
nodal  point  load.  It  too  resists  the  positive  displacement  of  node  3. 
Hence,  the  equation  of  motion  for  node  3  is 

-  m3d3  -  k3d3  -  Act  =  0 
or 

m3d3  +  k3d3  +  Acx  =  0 

Of  course,  node  2  does  not  have  the  shear  force,  so  its  equation  is 
m2d2  +  k2d2  =  0 

The  derivation  of  the  shear  stress,  x,  equation  can  be  seen  from 
Figure  10. 


Figure  10.  Diagram  for  Shear  Force  Equation  Derivation 


Fr  -  small  d3  compared  to  th.  the  small  angle  approximation  for  the 
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shear  strain,  Y,  applies,  i.e.  y  =  tan  y  =  d3/th,  and  x  =  G(t)  Y  = 
G(t)d3/th.  This  stress  is  applied  over  the  contact  area,  Ac. 
Therefore,  assembling  the  matrices  and  applying  the  boundary 
conditions  gives  the  elastic  system  equation  of  motion.  Taking  the 
Laplace  transform  of  the  system  and  applying  the  correspondence 
principle  to  the  shear  modulus  in  the  shear  pad  force  gives 


s2  M  D(s)  +  Ke  D(s)  =  -  G(s)K'v  D(s)  +U(s) 


(60) 


where 


D(s)  .  (  j 

\  D3(s)  j 

U(s)  -  /  0  ) 

( U(  s)  | 


pAL  4  1 

6  L 1  2 


=  2  -1 


L  L  -1  1  J 


K’  —  i  ®  ® 
r\  v  -  — 


th  Lo  1 


where  the  dj  values  are  displacements  of  the  ith  node.  The 
displacement  di  was  eliminated  by  the  zero  displacement  boundary 
condition.  M  is  the  assembled  mass  matrix,  Ke  is  the  assembled 
matrix  of  stiffness  corresponding  to  elastic  elements  (only  the  rods 
in  this  problem),  and  G(t)K'v  is  the  assembled  matrix  of  stiffness 
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due  to  viscoelastic  elements  (the  shear  pad).  Using  the  value  of 
G(s),  the  system  equation  of  motion  (EOM)  is 

[M  s^+Gi  K’vS1  2  +  (GoK'v  +  Ke)]  D(s)  =  U(s)  (EOM)  (61) 

Values  for  the  matrices  in  this  problem  are 

M  =  0.04660  4  1  kg 

.1  2  J 

Ke  =  6.894  x  10s  2  _1  N/m 

.  -1  1  . 

GoK'v  =  7600  0  0  N/  m 
.0  1. 

GiK’v  -  2950  0  0  N-sec1/2/m 

Lo  i . 


The  eigenvalue  problem  is  once  again  set  up  in  the  s1/2.piane 
by  defining  X  =  s1/2  ,  where  kjand  <t>j  are  the  jth  eigenvalue  and 
eigenvector.  This  gives 

M  +Gi  K'y  Xj  +  (Go  K'v  +  Kq)  ]  <>j( s)  =  0  (02’ 

which  can  be  solved  by  the  iteration  scheme  process  described  in 
Example  1.  The  first  four  eigenvalues  and  eigenvectors  are  shown  in 
Table  IV. 
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Table  IV.  First  Four  Open  Loop  Eigenvalue 
and  Eigenvectors,  Example  2 


Eigenvalue 

Eigenvector 

44.7541  ±  45. 4687i  , 

-44.7542  ±  44.0099i 

0.7162  ±  0.0088i 

1.0000 

[  0.6990  ±  0.0084i 

|  1.0000 

_ 

1  /2-Order  System  Formulation 

The  expansion  into  the  1/2-fractional  state  space  follows 
from  the  first  example: 


0  0  0  M 

s3  2  D(s) 

0  0  -M  0 

s3  2  D(s) 

s12 

0  0  M  0 

- 

s1  D(s) 

:  + 

0  -M  0  0 

s’  D(s) 

0  M  0  0 

s1  2  D(s) 

-M  0  0  0 

s1  2  D(s) 

_M  0  OKi. 

,  D(s) 

0  0  0  Kq  . 

,  D(s)  , 

0 
0 
0 
I  J 


U(S) 


(63) 


where 


K0  =  GqK'v  +  K 


e 


Ki  =GiK’v 
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Again  let  the  expanded  state  vector  be  denoted  by  an  underscore,  i.e. 
— * 

D(s).  Note  that  U(s)  is  a  scalar,  because  there  is  only  one  control 
force  in  this  example.  This  is  rearranged  to  the  appropriate 
s1  /2D(s)  =  A  D(s)  +  B  U(s)  form,  where 


0  0  0  M 

-i 

0  0  -M  0 

0  0  M  0 

0  -M  0  0 

0  M  0  0 

-MOOD 

M  0  0  Kt 

0  0  0  K0  . 

0  0  0  M 

0 

B  = 

0  0  M  0 

0 

0  M  0  0 

0 

M  6  0  Ki 

-7. 

(64) 


(65) 


U(s)  =  U(s) 

To  complete  the  system,  it  is  assumed  that  all  states  are  outputs 
and  that  there  is  no  feedforward  to  the  output,  i.e. 


Y(s)  =  CD(s)  +  DU(s) 


(66) 


where 


C  =  8  x  8  identity  matr  ix 
D  =  8  x  1  null  (zero)  matr  ix 


Open  Loop  Poles  Via  MATRIXv 

The  above  system  reduces  to  the  eigenvalue  problem  for  open 
loop  poles  when  the  control  vector,  U(s),  is  set  to  zero.  The 
eigenvalues  of  the  A  matrix  are  then  the  open  loop  poles.  The  first 
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eight  eigenvalues  are  found  by  using  the  EIG  function' on  A.  These 
are  shown  in  Table  V  with  the  closed  loop  poles.  Having  already 
found  the  firs:  four  by  the  iteration  solution  of  the  EOM,  the  second 
set  of  four  can  be  identified.  Notice  that  the  first  four  are  the  same 
as  found  by  the  iteration  method.  Unfortunately,  the  mode  shapes 
(eigenvectors)  for  the  second  set  must  be  found  by  modifying  the 
iterative  solution  of  the  EOM.  The  poles  are  shown  graphically  in 
Figure  11. 
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Figure  11.  Example  2  Open  Loop  Poles  for  Both  Vibrational  and 
Relaxation  Modes  1  and  2. 


Glased  Loop  Poles  Via  MATRIX^ 

The  system  is  posed  in  the  appropriate  form  for  the  Riccati 
equation  solver:  D(t)  =  A2  Q(t)  +  ABu(t)  .  There  are  two  weighting 
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matrices,  because  two  cases  are  studied  to  show  the  effect  of 
varying  weights.  This  is  also  required  to  get  a  significant  change  in 
the  pole  locations  -  an  example  of  the  usefulness  and  arbitrariness 
of  weighting  to  get  a  desired  pole  structure.  The  inputs  for  the 
Riccati  function  are 


A2  AB 
C  D 


NS  =  8 


MATRIXx  then  gives  KCi  for  RQi: 

KCi  = 

[-0.3709  -0.4680  14.986  1205.8  731.18  -3083.4  -3.1086D+06  1.5543D+06 

and  the  eigenvalues  of  the  system,  s1  QD{s)  =  (  A  -  BKCi  )  D(s)  ,  are 
found  by  applying  EIG(;  to  (  A  -  BKCi  ).  The  mode  to  which  each 
closed  loop  pole  corresponds  is  determined  by  marching  through 
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increasing  feedback  gain  and  observing  the  evolution  of  the  poles 
from  the  open  loop  ones.  The  poles  are  shown  in  Table  V. 

MATRIXx  gives  KQ>  for  RCfc: 

KC2  = 

[-0.3667  -0.4311  3376.1  33925  1  937  -3293  -2.6028D+0 8  1 .3014D+0 8  1 


The  corresponding  eigenvalues  are  also  shown  in  Table  V. 


Table  V.  s1/2-Plane  Open  and  Closed  Loop  Poles  for  Example  2 


HHSSIESBMHi 

Closed  Case  1 

Closed  Case  2 

44.7541  ±  45.4687i 

.  29.1442  ± 

59.8026i 

2.10  ± 

631. 92i 

44.7542  ±  44.0099i 

-26.8325  ± 

58.4049i 

-0.01  ± 

9.89i 

83.6001  ±  84.041  7i 

64.9435  ± 

95.061 8i 

54.06  ± 

77.67i 

-83.6000  ±  83.1  6 57 i 

-64.9542  ± 

93.3355i 

-54.08  ± 

77.71  i 

The  closed  loop  poles  are  shown  in  Figures  12  and  13. 
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Figure  12.  Example  2  Case  A  Closed  Loop  Poles  for  Both  Vibrational 
and  Relaxation  Modes  1  and  2. 
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Figure  13.  Example  2  Case  B  Closed  Loop  Poles  for  Both  Vibrational 
and  Relaxational  Modes  1  and  2. 


Comparison  of  Open  and  Closed  Loop  Poles 

The  stability  measures  based  on  the  s1/2-p!ane  open  and 
closed  loop  poles  are  shown  in  Table  VI,  and  the  oscillatory  motion 
parameters  are  shown  in  Table  VII. 
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Table  VI.  Stability  Measures  for  Example  2 


Mode  Re(O) 

Re(C) 

Re(0)-Re(C) 

Open  l/R 

Closed  l/R 

Case  1 

1  44.7541 

29.1442 

15.6099 

1 .0160 

2.0520 

-44.7542 

-26.8325 

-17.9217 

0.9834 

2.1766 

2  83.6001 

64.9435 

18.6566 

1 .0053 

1 .4638 

-83.6000 

-64.9542 

-18.6458 

0.9948 

1.4369 

Case  2 

1  44.7541 

2.1 

42.6541 

1 .0160 

300.91 

-44.7542 

-0.01 

-44.7442 

0.9834 

989.00 

2  83.6001 

54.06 

29.5401 

1 .0053 

1 .4367 

-83.6000 

-54.08 

-29.5200 

0.9948 

1  .4369 

Table  VII.  Example  2  Oscillatory  Motion  Parameters 


EV  sl/2  EV  si 

c 

con 

ati 

Open  Loop 

44.7541  ±  45.4687i  -64.473  ± 

83.6001  ±  84.041 7i  -74.031  ± 

4069. 8i 

14052. 7i 

0.0158 

0.0053 

4080.6 

14053.3 

4069 

14052 

Closed  Loop  Case  1 

29.1442  ±  59.8026i  -2727.0  ± 

64.9435  ±  95.061 8i  -481 9.1  ± 

3485. 8i 

12342. 3i 

0.6162 

0.3637 

4425.7 

13249.4 

3485 

1  2342 

Closed  Loop  Case  2 

2.1000  ±  631. 92i  -399318  ± 

54.0600  ±  77.67i  -3110.1  ± 

2654.1  i 

8397. 7i 

=1.0000 

0.3473 

399318 

8955.1 

2654 

8397 

Both  closed  loop  cases  show  improvement  in  stability  over  the  open 
loop  case.  Mode  one  damping  ratio  increases  from  0.0158  to  0.6162 
and  approximately  1.0.  Its  natural  frequency  also  increases  from 
4080.6  rad/sec  to  4425.7  rad/sec  and  399318  rad/sec,  in  cases  A 
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and  B  respectively.  Although  the  damping  ratio  increases  for  Mode  2, 
its  natural  frequency  decreases.  Interestingly,  the  damped 
frequency  for  both  modes  decreases.  Therefore,  resonance  occurs  at 
lower  frequencies  but  at  a  reduced  magnification  factor, 


12 


1  Q. 


| G(ico) |  =  -{ co/co n)  J  +  (2£co/con)  j  _  In  fact,  for  case  2  the  system 

is  approximately  critically  damped,  so  there  would  be  no  resonance 
for  that  mode.  Mode  2,  interestingly,  has  a  better  damped  response 
in  case  1  than  in  case  2  because  of  a  lower  G(ico)  and  faster  decrease 
in  the  amplitude  envelope.  At  resonance  the  larger  £  results  in  a 
lower  magnification  factor,  and  the  larger  £con  product  causes  the 
amplitude  to  die  faster  (e'£wnt)  .  Furthermore,  case  1  has  a  better 
separation  of  resonance  frequencies  than  case  2,  but  neither  has 
better  separation  than  the  open  loop.  Evidently,  minimizing  'ue 
performance  index  does  not  mean  that  all  stability  parameters  will 
improve.  It  is  possible  that  controlling  mode  1  is  having  an  effect 
on  mode  2  such  that  the  latter  is  being  destabilized.  Nevertheless, 
even  at  resonance  the  system  is  displaying  reduced  motion. 

Closed  loop  cases  1  and  2  demonstrate  how  weighting  can 
significantly  change  results.  In  the  modified  example  1  equal 
weighting  was  sufficient  to  create  a  significant  difference  in  the 
pole  structure,  while  here  in  example  2  weights  of  106  and  109 
were  reauired  to  make  significant  changes.  Hence,  the  guideline  for 
weighting  is  the  desired  pole  structure  at  the  end.  Furthermore, 
various  combinations  of  weights  can  result  in  similar  pole 
structures.  For  example,  weighting  velocity  at  106  and  control  at 
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10-3  yields  results  similar  to  weighting  velocity  at  109.  Fora 
guideline  on  weighting  remember  from  example  1  that  the  velocity 
and  derivative  of  order  3/2  in  feedback  have  strong  effects  on  poles, 
while  the  derivatives  of  lesser  order  have  much  less  effect. 
Therefore,  to  make  a  significant  change  in  pole  location,  one  should 
weight  the  velocity  or  derivative  of  order  3/2  ;  while  weighting 
position  should  only  make  small  changes  to  the  pole  structure.  Since 
the  control  vector  becomes  full  state  feedback,  weighting  it  has 
mixed  effects,  because  it  includes  all  derivatives.  However, 
example  2  results  indicate  that  the  conrol  vector  effect  is 
dominated  by  the  higher  order  derivatives,  which  allowed  weighting 
velocity  at  106  and  control  at  10-3  to  yield  results  similar  to 
weighting  velocity  at  109- 
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V. 


Conclusions 

Quadratic  control  theory  can  be  used  to  find  an  optimal 
feedback  control  law  at  steady  state  conditions  for  a  structure 
incorporating  a  fractional  calculus  model  of  viscoelasticity. 

Following  linear  regulator  control  theory,  it  has  been  shown  that 
under  large  time  conditions  the  optimal  control  law  is  linear  state 
feedback,  assuming  an  optimal  time-varying  gain  matrix,  K,  which  is 
asymptotically  constant  for  large  time,  for  the  minimized 
periormance  index.  A  Riccati  equation  is  created,  such  that  a 
fractional  derivative  system  of  order  1/n  (n  an  integer)  can  be 
represented  by  an  equivalent  first-order  system  in  Riccati  equation 
solving  routines.  Therefore,  current  tools  can  be  used  to  easily 
solve  problems  of  this  kind.  By  expanding  the  structure's  equation  of 
motion  into  a  fractional  state-space  system,  the  control  theory  can 
be  applied  to  structures  incorporating  both  passive  and  active 
damping.  This  has  been  illustrated  in  two  example  problems. 

The  theory  is  limited  to  optimal  gain  matrices  which  are 
asymptotically  constant  for  large  time.  The  optimal  control  ,  as 
determined  by  linear  regulator  theory,  is  asymptotically  linear 
feedback  as  time  tends  to  infinity.  The  evaluation  of  the  Hamilton- 
Jacobi-Bellman  equation  cannot  produce  a  Riccati  equation  in 
general  time  for  a  time-varying  optimal  gain  because  of  coupling  of 
the  gain  matrix  and  the  state  vector  due  to  fractional  derivatives 
acting  on  the  control  vector.  Only  when  time  tends  to  infinity  can 
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these  derivatives  be  neglected.  Then  the  optimal  gain  matrix  can  be 
decoupled  from  the  state  vector,  and  an  algebraic  Riccati  equation 
can  be  identified  for  time  large  and  for  an  optimal  gain  matrix  which 
is  asymptotically  constant  for  large  time. 


I 
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VI.  Recommendation  for  Further  Work 

The  next  step  is  to  test  the  theory  by  experiment.  Captain 
Kevin  Klonoski's  thesis  research  (16)  has  developed  a  prototype 
sensor  for  fractional  derivative  motion,  and  a  design  useful  for 
experimentation  should  be  developed.  Response  of  a  simple  beam 
structure  should  then  be  predicted  using  the  theory  and  compared  to 
experimental  results  using  such  a  sensor  to  feed  back  the  fractional 
state  vectors.  This  will  also  require  the  development  of  a  discrete¬ 
time  model  and  theory,  if  digital  control  is  to  be  used.  Thus, 
verification  requires: 

1.  Design  of  a  fractional  derivative  sensor  suitable  for 
experiment  instrumentation; 

2.  Development  of  a  discrete-time  version  of  the  theory  in 
order  to  apply  digital  control  (if  desired); 

3.  Comparison  of  the  predictions  with  experimental  data  for  a 
simple  system. 

Further  theoretical  research  could  examine  the  effect  of 
optimal  performance  index  forms  other  than  that  used  in  linear 
regulatory  theory.  There  is  evidence  that  time  delay  in  feedback 
loops  can  cause  exponential  instability.  Hence,  further  research 
could  siudy  the  effect  of  time-delayed  feedback  on  systems  using 
this  theory. 

Another  possible  area  of  interest  is  the  investigation  of  the 
effect  of  this  theory  on  modes  over  a  large  range  of  frequencies.  I 
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have  only  worked  with  the  first  and  second  modes  and  there  are 
significant  differences  in  results  between  the  two  example 
problems.  Examination  of  the  effect  on  a  larger  range  of  modes  and 
comparison  with  physical  theory  on  how  those  higher  modes  should 
be  effected  could  give  insight  to  the  accuracy  and  utility  of  this 
control  model. 
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Appendix  A:  Reaularlv-Varvina  Functions 

The  presentation  in  this  appendix  is  excerpted  with  minor 
changes  from  Warhola  (26). 

"Regularly-varying"  functions  are  those  functions,  f(t),  for 
which  f(ct)/f(t)  approaches  a  finite  non-zero  limit  as  t  ->  00  ;  for 

any  such  function,  the  limit  is  necessarily  of  the  form  cP  as  a 
function  of  the  parameter  c,  with  .  00  <  p  <  °o  (10).  In  this  limit,  f 

has  the  asymptotic  behavior: 

f(ct)A^  cp  f(t),  t  -»  00  (A.1) 

For  such  a  function,  the  notation  P(f)  =  p  indicates  that  p  is  the 
"power"  of  f. 

If  Eq  (A.1 )  holds  with  p  »  0,  the  function  is  "slowly-varying." 
To  suggest  a  logarithm,  L(t)  is  used  to  denote  any  otherwise 
unspecified  slowly-varying  function.  Then 

L(ct)A/  L(t)  ,  t->°°  (A. 2) 

Any  regularly-varying  function  f(t)  can  be  expressed  in  the  form 

f(t)  =  L(t)  tP  (A. 3) 

For  many  purposes,  L  can  be  treated  as  if  it  were  a  constant,  even 
though  it  may  converge  to  zero  or  diverge  to  infinity.  For  example, 
powers  and  products  of  regularly-varying  functions  are  also 
regularly-varying,  with  the  obvious  exponents.  In  particular,  any 
function  which  approaches  a  non-zero  constant  value  is  slowly- 
varying. 
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If  the  slope  on  a  log-log  graph  of  a  differentiable  function, 
f(t),  approaches  a  constant,  p,  for  t  -»  00 , 

d(!n  f(t))/d(ln  t)  =  tf'(t)/f(t)  p,  t  ->  oo  (A. 4) 

then  f  is  regularly-varying  with  P(f)  =  p.  Additionally,  f  is  also 
regularly-varying  (unless  f  is  merely  slowly-varying),  with 
P(f’)  =  p-1 : 

f  (t)  p  f(t)/t  ,  p  *  0  ,  t  oo  (A. 5a) 

f  (t)  =  o(  f(t)/t  ),  p  =  0  ,  t  ->  oo  (A. 5b) 

The  exception  for  p  =  0  occurs  because  when  f  is  approaching  a 
constant,  f'  may  be  approaching  zero  much  faster  than  1/t. 

On  the  contrary,  when  f  is  regularly-varying,  its  indefinite 
integral  F  is  always  regularly-varying  (10).  When  F  diverges  as 

t  — >  oo  , 


F(t)  =  f  f(x)  dx~  t— ,p  =  P(f)  >  -1  ,  t — ►oo 

Jo  P+1  (A. 6) 

just  as  if  f  were  actually  a  power.  When  p<-1  ,  the  integral  in 
Eq  (A. 6)  approaches  a  constant  and  the  tail  of  the  integral  is 
regularly-varying: 

f  f  (x)  dx  ~  ,  p  =  P(r)  <  -1  ,  t  ->  oo 

it  IP+1I  (A. 7) 


The  result  analogous  to  Eq  (A5.b)  is 
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t  f(  t)  «  o(F)  ,  p  =  P(f)  <  -1  ,  t->oo 


A  useful  result  is  added  here  for  the  convolution  of  two 


(A. 8) 


regularly-varying  functions  of  a  certain  class.  If  f(t)  and  g(t)  vanish 
for  t  <  0,  and  are  regulariy-varying  at  infinity  with  powers 
P(f)  =  p  >  -1  and  P(g)  =  q  >  -1 ,  then 


r\ 

f(  t-T)  g(-c)  dr—  t  f(  t)g(t)  ,  p,  q 

lo  (P-«I+1)! 


>  -1  ,  t  -»  OO 


(A. 9) 


where 


tD  e_t  dt  ,  p  >  -1 
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This  is  a  listing  of  the  MATRIXx  command  file  (14:7-3  -  7-6) 
called  4BEAMEVS.EXC.  This  file  determines  the  open  loop 
eigenvalues  and  eigenvectors  for  the  first  mode  of  a  fourth-order 
beam  equation  of  motion  which  incorporates  a  fractional  calculus 
model  of  viscoelasticity  of  order  1/2,  i.e.  the  fractional  derivatives 
are  of  order  1/2. 

The  file  is  interactive.  The  user  must  supply  the  assembled 
mass,  elastic  stiffness  (Ko),  and  viscoelastic  stiffness  (Ki) 
matrices,  as  defined  in  Chapter  IV.  The  user  must  also  supply  an 
initial  guess  of  the  eigenvalue  and  eigenvector.  The  iteration 
process  is  not  automated,  so  the  file  will  as!-'  which  of  four  possible 
roots  should  it  follow.  Choose  one  set  of  signs  to  follow  at  a  time 
until  the  root  converges  within  an  acceptable  error.  Complete  this 
four  times  to  get  all  four  eigenvalues  and  eigenvectors.  Answer  the 
query  on  which  root  to  follow  by  RT(element  number  of  proper  root). 
The  vector  of  roots  is  called  RT,  so  RT(#)  specifies  a  specific 
element. 

The  eigenvalues  and  eigenvectors  are  written  to  file  EV.DAT  as 
vector  EVAL  and  matrix  EVEC.  Each  column  of  EVEC  is  an 
eigenvector. 

DISPLAYOEIGENVALUE/VECTOR  SOLN  FOR  4TH  ORDER  BEAM  THEORY) 
INQUIRE  N  ’ENTER  NUMBER  OF  DOF 

INQUIRE  M  ’ENTER  ASSEMBLED  CONSISTENT  MASS  MATRIX’ 

INQUIRE  KO  ’ENTFR  ASSEMBLED  KO  MATRIX’ 

INQUIRE  KI  ’ENTER  ASSEMBLED  KI  MATRIX’ 


El  K=K1 ; 

INQUIRE  LI  'ENTER  INITIAL  UVMBDA1  GUESS’ 

INQUIRE  EV1  'ENTER  INITIAL  EIGENVECT0R1  GUESS’ 

K2=L1*E1  K+EOK; 

K2INV=INV(K2); 

EV1  ERR=1; 

WHILE  EV1cRR>5D-6,... 

LHS=K2INV*M*EV1 
LHSMAX=LHS(1 );... 

FOR  l=2:N,  IF  ABS(LHS(I))>ABS(LHSMAX),  LHSMAX=LHS(I);END,END,.. 
EVT=LHS/LHSMAX,... 

EV1  ERR=ABS(EVT(1)-EV1  (1));... 

FOR  U2:N,... 

IF  EVT(l)oEV1  (I),... 

ERR=ABS(EVT(I)-EV1  (I)),... 

IF  ERR>EV1  ERR.EV1  ERR=»ERR;END,... 

END,... 

END,... 

IF  EV1  ERR>5D-6,  EV1=EVT;END,... 

EM) 

LHSMAX=LHSMAX 

H=1/LHSMAX; 

P=[1 ,0,0,0  ,H]; 

RTBASE=ROOTS(P) 

FOR  J=1 :4,... 

L1=RTBASE(J);... 

CHOICER;... 

OLDRLERR=1 ;... 

OLDIMERR=1;... 

ITER=1  ,... 

WHILE  CHOICER,... 

K2=L1*E1K+E0K;... 

K2INV=INV(K2);... 

EV1ERR=1 ;... 

WHILE  EV1  ERR>5D-6,... 

LHS=K2INV*M*EV1 ;... 

LHSMAX=LHS(1 );... 

FOR  l=2:N,  IF  ABS(LHS(I))>ABS(LHSMAX), 
LHSMAX=LHS(I);END,END,... 

EVT=LHS/LHSMAX;... 

EV1  ERR=ABS(EVT(1  )-EV1  (1 ));... 


FOR  l=2:N,... 

IF  EVT(l)oEV1  (I),... 

ERR=ABS(EVT(I)-EV1(I));... 

IF  ERR>EV1  ERR.EV1  ERR=ERR;END,... 

END,... 

END,... 

IF  EV1ERR>5D-6,  EV1=EVT;END,... 

END,... 

LHSMAX=LHSMAX,... 

H=1/LHSMAX;... 

P=[1 ,0,0,0,  H];.. . 

RT=ROOTS(P),... 

INQUIRE  LI  TEMP  'WHICH  ROOT  TO  CONTINUE  ON  BRANCH?',... 
ERRTEST=L1  TEMP-L1  ;... 

NEWRLERR=REAL(ERRTEST);... 

NEWIMERR=IMAG(ERRTEST);... 

DRLERR=NEWRLERR-OLDRLERR,... 

DIMERR=NEWIMERR-OLDIMERR,... 

INQUIRE  CHOICE  'CONTINUED  STOP=1’,... 
OLDRLERR=NEWRLERR;... 

OLDIMERR=NEWIMERR;... 

L1=L1TEMP,... 

ITER=ITER+1,... 

END,... 

EVAL(J)=L1TEMP,... 

FOR  1=1  :N,  EVEC(I,J)=EV1  (l),END,... 

BVD 

PRINTCEV.DAT’,  EVAL) 

P  R I  NT(’  EV .  DAT ,  E  VEC) 

RETURN 
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The  approach  requires  expanding  the  structure's  equation  of  motion  into  a  fractional 
state-space  system  of  order  1/n,  where  n  is  an  integer  based  on  the  viscoelastic 
damping  material  constitutive  law.  This  approach  restricts  the  theory  to  materials 
which  have  rational,  fractional  derivatives  in  their  constitutive  laws.  An 
equivalent  first-order  system  is  then  formed  and  used  to  derive  the  optimal  control 
theory  in  linear  regulator  problems.  The  quadratic  performance  index  used  in  linear 
regulator  theory  is  used,  and  equations  similar  to  those  derived  for  linear  regulator 
problems  using  first-order  systems  are  developed.  The  optimal  control  law  is  asymp¬ 
totically  linear  feedback  of  the  state  vector  for  time  large.  An  equation  which 
defines  the  optimal  gain  matrix  for  the  feedback  is  derived  and  is  asymptotically 
an  algebraic  Riccati  equation  for  long  time  and  for  gain  matrices  which  are 
asymptotically  constant  for  large  time.  Since  an  algebraic  Riccati  equation  can  be 
defined,  current  solving  routines  can  determine  the  asymptotically  constant  optimal 
gain  matrix  for  large  time  for  the  fractional  state-space  system. 

In  the  general  time  case,  no  Riccati  equation  can  be  derived  because  of  coupling 
between  the  optimal  gain  matrix  and  the  state  vector  due  to  fractional  derivatives 
of  the  control  vector.  Only  when  gain  matrices  are  asymptotically  constant  for 
time  large  and  time  is  large,  do  they  uncouple. 

The  theory  is  illustrated  by  two  examples.  A  simply-supported  viscoelastic  beam 
with  controllers  illustrates  the  solution  process  in  Example  1.  The  beam  example 
incorporates  the  fractional  calculus  viscoelastic  behavior  in  a  structure  element, 
while  an  axially  deforming  rod  in  Example  2  incorporates  the  behavior  through 
a  viscoelastic  shear  force  applied  at  a  node  by  a  damping  pad.  The  example 
equations  of  motion  are  numerically  solved  using  the  commerc ially-ava ilable 
control  analysis  software  package  called  MATRIX^,. 


